NASA CONTRACTOR 
REPORT 



CM 





nm&efm s 



p^rup?* 



LOAN COPv. 
A^-WL TECHKJCAL utV-° 

K'RTLAND AF^^r 






DESAP 1 - A STRUCTURAL DESIGN PROGRAM 
WITH STRESS AND DISPLACEMENT CONSTRAINTS 

Volume I: Theoretical and User's Manual 



/. Kiusalaas and G. B. Reddy 

Prepared by 

THE PENNSYLVANIA STATE UNIVERSITY 

University Park, Pa. 16802 

for George C. Marshall Space Flight Center 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • FEBRUARY 1977 



TECH LIBRARY KAFB, NM 

IMIHII 



1. REPORT NO. 

NASA CR-2794 



2. GOVERNMENT ACCESSION NO. 



3. RECIPIENT'S u^ 



QOblBS? 



TITLE AND SUBTITLE 

DESAP 1 - A Structural Design Program with Stress and Dis- 
placement Constraints. Vol. I: Theoretical and User's Manual 



S. REPORT DATE 

February 1977 



6. PERFORMING ORGANIZATION CODE 



7. AUTHOR(S) 

J. Kiusalaas and G. B. Reddy 



8. PERFORMING ORGANIZATION REPORT # 
M-202 



PERFORMING ORGANIZATION NAME AND ADDRESS 



10. WORK UNIT NO. 



Department of Engineering Science and Mechanics 
The Pennsylvania State University 
University Park, Pennsylvania 16802 



1 1. CONTRACT OR GRANT NO. 

NCA-8-71/72 



12. SPONSORING AGENCY NAME AND ADDRESS 



National Aeronautics and Space Administration 
Washington, D. C. 20546 



13. TYPE OF REPORV & PERIOD COVERED 

Contractor Report 



14. SPONSORING AGENCY CODE 



IS. SUPPLEMENTARY NOTES 



This work was done under the sponsorship of the Structures and Propulsion Laboratory 
of the Marshall Space Flight Center, Alabama. 



16. ABSTRACT 

DESAP 1 is a finite element program for computer-automated, minimum weight design 
of elastic structures with constraints on stresses (including local instability criteria) and 
displacements. Volume 1 of the report contains the theoretical and user's manual of the pro- 
gram. Sample problems and the listing of the program are included in Volumes 2 and 3, 
respectively. 

The static analysis portion of DESAP 1 is based on the SOLID SAP finite element pro- 
gram developed at the University of California, Berkeley. In design, the stress ratio method 
is employed for the stress constraints, whereas the displacement constraints are handled by 
solving the appropriate optimality criterion. 

The element subroutines have been organized so as to facilitate additions and changes 
by the user. As a result, a relatively minor programming effort would be required to make 
DESAP 1 into a special -purpose program to handle the user's specific design requirements 
and failure criteria. 

DESAP i is a companion program of DESAP 2, "A Structural Design Program with 
Stress and Buckling Constraints. " With the exception of a few cards, the same data deck can 
be used for both programs. 

This is Volume 1 of three volumes. 



17. KEV WORDS 



18. DISTRIBUTION STATEMENT 



STAR Category 39 



19. SECURITY CLASSIF. (of thl» r«pcrtl 

Unclassified 



20. SECURITY CLASSIF. (of thl. page) 

Unclassified 



21. NO. OF PAGES 

122 



22. PRICE 



$5.50 



*For sale by the National Technical Information Service, Springfield, Virginia 22161 



ACKNOWLEDGEMENTS 

The development of the program was sponsored by the National 
Aeronautics and Space Administration, George C. Marshall Space Flight 
Center, Huntsville, Alabama, under cooperative Agreement NCA-8-71/72 
with the Pennsylvania State University, University Park, Pennsylvania. 
The technical monitor of the work was J. E. Key. 

The authors would also like to acknowledge the contributions of 
M. Isreb, who assisted in the development of the shear panel element 
and the stress-constrained redesign subroutines. 



11 



TABLE OF CONTENTS 
VOLUME 1 
Abstract 
Acknowledgements 

A. INTRODUCTION 

B. THEORETICAL BACKGROUND 
B.l. Element Properties 
B.2. Stress Constraints 

B.3. Displacement Constraints 

B.4. Stress and Displacement Constraints 

B.5. Uniform Scaling Operation 

C. ORGANIZATION OF THE PROGRAM 
C.l. Overlay Tree 

C.2. Element Subroutines 

C.3. Storage Requirements 

C.4. Requirements for Auxiliary Storage Devices 

C.5. Assembly and Solution of Equations 

C.6. Cutoff Criteria and Restart Option 

D. DESCRIPTION OF INPUT DATA 
D.l. General Input Data 

E. BAR ELEMENT 

E.l. General Information 
E.2. Construction Code No. 1 
E.3. Construction Code No. 2 

F. BEAM ELEMENT 

F.l. General Information 
F.2. Construction Code No. 1 
F.3. Construction Code No. 2 

G. PLANE STRESS QUADRILATERAL OR TRIANGULAR ELEMENT 
G.l. General Information 

G.2. Construction Code No. 1 
G.3. Construction Code No. 2 

H. QUADRILATERAL SHEAR PANEL 
H.l. General Information 
H.2. Construction Code No. 1 

I. PLATE QUADRILATERAL OR TRIANGULAR ELEMENT 

1 . 1 . General Information 

1.2. Construction Code No. 1 

J. BOUNDARY ELEMENT 

References 

iii 



A.l 
A. INTRODUCTION 

DESAP I is a finite element program for automated design 
(synthesis) of linear-elastic structures under static loads. The 
design objective is to find the element sizes (cross-sectional areas, 
plate thicknesses, etc.) that minimize the total structural weight. 
The layout of the structure is not changed during the design procedure. 

The primary constraints used in the synthesis algorithm are 
upper limits on stresses and displacements. The stress limits may be 
prescribed in the form of yield criteria, local instability criteria, 
or both. 

The program also allows the use of secondary constraints , which 
consist of minimum allowable element sizes and size proportion con- 
straints (the requirement that the sizes of specified elements be 
equal, or have a prescribed ratio). 

The design procedure is based on the assumption that the loading 
is independent of the element sizes (dead loading of prescribed magni- 
tude) . Allowance has been made for size-dependent loads, e.g., thermal 
stresses and gravity loading, but these loads must be used with special 
precautions . 

An iterative process is used as the synthesis algorithm, each 
iterative step consisting of three parts: 

1) Analysis of the current design. 

2) Evaluation of the analysis. 

3) Redesign of the structure, based on the results of 2) and 3). 
The analysis portion of the design procedure is essentially the 

SOLID SAP finite element program developed by E. L. Wilsonfl]. It 
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was necessary to make numerous modifications to the original program, 
most of which were dictated by the special requirements of the redesign 
operations. Other major changes are: 

1) A shear panel quadrilateral element was added. 

2) Some element subroutines were reorganized in order to increase 
computer efficiency and eliminate unnecessary portions of the 
program. 

3) Temperature-dependent material properties were incorporated for 
all elements, except the beam. 

4) Certain parts of the program were converted to double precision 
computation. 

The classical stress ratio method is employed in the redesign with 
respect to stress constraints . This procedure will drive the final 
design to the fully stressed design, which does not necessarily coincide 
with the minimum weight distribution of the material. Unfortunately, 
at the present state-of-the-art, more rigorous design methods are en- 
tirely impractical for structures of reasonable size, since they re- 
quire the inversion of the structural stiffness matrix, or the use of 
numerical search techniques . 

For the displacement constraints , the redesign procedure described 
in [2] is used. The redesign formula is obtained from the minimum 
weight criterion; consequently, the element sizes are driven towards 
an optimal weight design. 

As pointed out repeatedly in existing literature, neither the fully 
stressed design, nor the optimal weight design are necessarily unique, 
i.e., they have a local rather than a global character. It follows 
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therefore, that the choice of the initial design plays an important role 
in determining the design to which the synthesis algorithm converges. 
Numerous examples seem to indicate, however, that the weight differences 
between the various converged designs are slight, although there may be 
large differences in the distribution of the material. 

In the presence of both stress and displacement constraints, each 
iterative redesign operation first determines the element sizes from 
the stress ratio method. The results are then used as minimum size 
constraints during the application of the displacement-constrained re- 
design formulas. 

The basic structural elements that have been presently programmed 
are: 

1) Three-dimensional bar. 

2) Three-dimensional beam. 

3) Quadrilateral shear panel. 

4) Anisotropic, plane stress quadrilateral and triangle. 

5) Anisotropic, quadrilateral and triangular plate. 

6) Boundary elements. These are not subjected to redesign, but 
are used to represent elastic supports, or to apply constraints. 

A very important feature of the program is the organization of 
the element subroutines such that they can easily be adapted to the 
user's special requirements with the smallest possible programming 
changes. In making this provision, we recognized the fact that it is 
virtually impossible to create an all-purpose synthesis program. The 
reason is that structural design requires a much larger and more varied 
volume of input information than analysis. In the case of beam elements, 
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for example, the cross-sectional area, the two moments of inertia and 
the torsional constant suffice for the purposes of analysis. In de- 
sign we must add the section moduli, buckling information and the 
yield criterion; in addition we must specify how all these properties 
vary with the design variable. Since each design situation may use 
elements of special construction and shape, and different design cri- 
teria, it is clearly impractical to make provisions for all the possi- 
ble forms of design data in a single program. 

Finally, it must be pointed out that DESAP 1, despite its flexi- 
bility, can still handle only a limited amount of design information. 
Consequently, the results of the program are to be taken as a pre- 
liminary design in the sense that it gives the desired proportions of 
the structure, but the structural details must still be determined by 
the designer through conventional design practices. 
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B. THEORETICAL BACKGROUND 



B. 1 Element Properties 

th 
The size of a typical (i ) element is denoted by A. . It repre- 
sents the cross-sectional area or the panel thickness for one and two- 
dimensional elements, respectively. The weight of the element can be 
written as 



W i = p i A i ' i = 1, 2, .... I (B.l.l) 

where p. is the unit weight of the element. If the minimum cost, rather 
than weight is the design objective, then p. should be taken as the 
unit cost. 

It is seldom desirable to have each A. as an independent 
design variable. Equal size constraints can be imposed by introducing 
the design variables D, m=l, 2 , ...,M, where M _< I , which are in- 
dependently variable, and expressing each element size in the form 



A. = n-D , (B.1.2) 

l l m J 

where r\. is called the design variable fraction of the element. Both 
r). and m must be specified for each element. Equal sizing of two 
elements is obtained, for example, by prescribing the same m for each 
element together with r|. = 1. In addition, the scheme permits propor- 
tional sizing if the same m, but different values of n- are used. 

The stiffness matrix of each element is restricted to the form 

[K,] =[k. (1) ]A. + [k. (2) ]A/ CB.1.3) 



B.1.2 

where [K.] is the element stiffness matrix, [k. J ] and [k. J ] are the 
unit stiffness matrices of the element, and n . is the inertia exponent . 
The first part of (B.1.3) represents the action of direct stresses, 
whereas the second part is due to bending or torsion. The value of n . 
depends on the physical nature of the design variable. For example, 
n. = 1 for thin-walled beams if the wall thickness only is being varied. 
If all dimensions of the cross section are scaled uniformly, 
then n. = 2. For plates, where the thickness is subject to design, 

we have n. =3. The unit stiffness matrices of each element are stored 

i 

separately on an external storage device, together with m, r\. and n . . 

It should be noted that n. i s determined by the relationship between 

the size and the moment of the inertia of the element: 

n. 

I. = j.A. 1 , (B.1.4) 

l J i l ' v ' 

where i . is the unit moment of inertia. 

J i 

The vector of 'internal forces {N. } of an element is recovered from the 

l 

element nodal displacements {u. } by the formula 

{N ± } = [S i ]{u i > + {T ± } . (B.1.5) 

If the element stiffness matrix is given by (B.1.3), then the force 
recovery matrix [S.] has the same form: 

[S.] = [s[ l) ]\ + [sf V. 1 , (B.1.6) 

where [s. ] and [s. ] are the unit force recovery matrices . The 
force vector {T. } also consists of two parts: 

{f.} = {t. (1) }A i + {t. (0) } . (B.1.7) 
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{t. J } being the unit thermal force vector and {t. J } represents the contri- 
bution of size-independent loading; e.g., fixed-end forces of beam ele- 
ments due to dead loading. Again, [s. ], [s. ], {t. } and {t. } 
are stored for each element. 

The element load vector {Q. } is similarly separated into two components: 

{Q.} = {q i C1) }A i + {q. C °>} , (B.1.8) 

where {q. } is the unit load vector due to size -dependent loading (thermal 
and gravity loads), and {q. } represents the contribution of dead loads; 
both vectors are stored for each element. 

The form of equations (B.1.7) and (B.1.8) allows only for a uniform 
temperature increase within an element, i.e., it assumes that thermal 
expansion causes extension without bending. This excludes, for example, 

the effects of thermal gradients through the thickness of plates, for which we 

j-2-v n - C2i n i 
would require additional terms of the type {t. }A. and {q. }A. 



for {T. } and {Q. }, respectively. 
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B.2 Stress Constraints 

At the present time, the only practical means of handling stress 
constraints for large structures seems to be the concept of fully 
stressed design (FSD) . In a fully stressed structure each element 
reaches its maximum permissible stress level under at least one load 
condition, unless it is governed by the minimum size constraint. FSD 
does not generally coincide with the minimum weight design, except for 
statically determinate structures, but the differences in weight are 
small in most cases. 

A typical failure criterion of an element can be expressed in the 
general form 

f({N i >, {N*},A.) = 1 , (B.2.1) 

* 
where {N. }is the internal force vector of the element and {N.} contains 

the allowable forces. Equation (B.2.1) may represent a criterion for 

any kind of failure, such as yielding, fracture or local instability. 

Let A. be the element size for the current design, and A. the size 
1 b 1 

predicted for next (improved) design. The corresponding internal forces 
are denoted by {N. } and {N.}, respectively. For the sake of clarity, 
we limit the discussion at this time to a single load condition, and 
assume that each element is subject to a single failure criterion. 

If the improved design is to be fully stressed, it must satisfy 

f({N!}, {N?},Ap = 1 . (B.2. 2) 

The main difficulty in using (B.2. 2) is that it requires a knowledge of 
{N.}, i.e. the changes in the nodal forces of each element caused by the 
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redesign. This information can be acquired only by an inversion of the 
structural stiffness matrix, which is impractical for large structures, 
since the banded form of the structural stiffness matrix is lost upon 
its inversion. 

The problem is greatly simplified if we assume that the nodal forces 
do not change upon redesign. Equation (B.2.2) then becomes 

f({N.}, {N*},a|) = 1 , (B.2.3) 

t 

which can immediately be solved for A. , one element at a time. We have 

1 

now arrived at the stress ratio method of redesign, which is adopted 
for DESAP 1. 

The assumption of no change in the nodal forces is valid only for 
statically determinate structures under dead (size-independent) loading, 
in which case a single redesign will result in FSD. For all other pro- 
blems, equation (B.2.3) is an approximation, and it must be applied 
iteratively, updating {N.} each time, before FSD is reached. If size- 
dependent loading is present (e.g., gravity or thermal loads), the 
approximation may be poor and cause difficulties of convergence. 

When several load conditions and design criteria are used, the 
element size is calculated for each combination of loading and design 
formula, and the largest value is chosen as A. . 

Once the new size of the element is known, the corresponding de- 

i i 

sign variable is determined by (B. 1.2): D = A./n., or the minimum 
& m 1 ' '!> 

i * * 

size constraint: D = D (D is the prescribed lower limit), whichever 

m m m c J ' 

is larger. If the design variable D is common to several elements, 
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i.e., if equal or proportional size constraints are used, the maximum 

i 

value of D is used: that is 
m ' 

it i * 

D = max A./n. or D = D , fB.2.4) 

m . 1 1 mm' K J 

lem 

whichever is greater. The notation iem indicates that maximum is ob- 
tained from all elements i that share the design variable number m. 
The ratio 

R = D /D 
m mm 

is called the stress ratio of the design variable. In order to 
monitor the progress of the design sequence, the maximum and minimum 
stress ratios of the structure, 

R = max R and R . = min R , (B.2.6) 

max m mm mm v J 

m 

i 
are computed after each redesign cycle. The new design {D } is said to 

be a stress -critical design if 

1 - 5 < R <l+6, (B.2.7) 

— max — ' *• J 

where 6 is a small parameter prescribed by the user. Similarily the 
design is considered to be fully stressed if, in addition to being 
critical, it also satisfies 



1 - 6 < R . < 1 + 6 . (B.2.8) 

-■ mm — J 

The values of R , R . and the corresponding design variable and 
max mm r & & 

load numbers are printed out after each redesign cycle. 
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B.3 Displacement Constraints 

Displacement constraints are handled in DESAP 1 by the technique 
described in Reference [2] . The procedure is not only applicable to 
displacements, but can also be used for buckling or frequency constraints 
see, for example, Reference [3]. 

Constraints on nodal displacements (or rotations) have the form 

u r < u* , r = 1, 2, ..., R , (B.3.1) 

* 
where u is a nodal displacement, u its maximum allowable value, and R 
r r 

denotes the number of constraints. The constraints can be divided into 
two categories: if the equality sign governs the optimal design, the 
constraint is said to be active ; if the inequality occurs, the con- 
straint is passive . This division is not known, of course, until the 
final design is reached, but for notational convenience we presume that 
the active constraints are listed first, so that (B.3.1) can be re- 
placed by 



* 



u r = u r , r = 1, 2, ..., R act , 



(B.3. 2) 



* 



u r <U r ' r = R act + 1 ' "•' R '' 
where R is the number of active displacement constraints. 

aCt 

In addition to displacement constraints, we must also account for 
limits on element sizes: 

A. > A. , i.e. , D > D . (B.3. 3) 

l—i ' m — m *■ ■* 

The design variables are also deemed as either active or passive, de- 
pending on their role in satisfying the displacement constraints. 
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A passive design variable is governed by the minimum size constraint 

* 

in the final design, i.e., D = D , whereas an active design variable 
6 ' m m - 

* 

is determined by the displacement constraint, in which case D' > D . 

3 r mm 

In the mathematical treatment that follows, it is convenient to 
replace (B.3.1) and (B.3.3) by the equality constraints 

u + a 2 = u* , D - b 2 = D* , (B.3.4) 

rrr m m m 

where a and b are to be viewed as variables free from constraints, 
r m 

The design objective is to minimize the total weight of the struc- 
ture W = EW. . Substituting for W. from (B.l.l) and (B.1.2), we have 
i 1 1 

W = Ep.n-D . (B.3.5J 

Minimizing (B.3.5J subject to (B.3.4) is equivalent to making the 
following function stationary: 

V= ?p n.D + EA (u + a 2 ) - Eu m (D -b 2 ) , (B.3.6) 

where A and y are non-negative Lagrangian multipliers . The operations 

3V/3D = 0, 3V/Sa = and 3V/3b = yield, respectively, 
m r m ' r J ' 

E p.n- +EAu - y D =0, (B.3.7) 

iem x i r r r ' m m m 

A r a r = , (B.3.8J 

y b = , (B.3.91 

mm J 

where we used the notation ( ) m = 3( )/3D . 

,m m 
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Inspection of (B.3.4), (B.3.8) and (B.3.9) reveals that 



A 
r 



% 



* 

> if u = u [active constraints) 
— r r 



* 



=0 if u < u (passive constraints) , 



* 

> if D = D (passive design variables) 

— m m r 6 J 

* 

= if D > D (active design variables) , 
mm 6 



which enables us to rewrite the optimality criterion (B.3.7) as 

f ■> if D = D* 

Z p.n- + E A u -i m * (B.3.10) 

iem X 1 r act r r,m = if D > D 

L mm 

The notation r act shows that the sum is to be taken over the active 

displacement constraints. Introducing the unit weight of the m 

design variable 

v (B.3.11) 

p m = Y p i n i • 
ism 

(B.3.10) becomes for the active design variables 

1 R act 

— Z A u = -1 (B.3.12) 

p . r r.m v * ' 

m r= 1 

The redesign formula for the active design variables is obtained 

directly from the optimality criterion. Multiplying both sides of 

(B.3.12) by (l-a)D , where a is a constant to be determined later, and 

rearranging terms , we get 



D R act 

D = oD - (1-cO— Z A „ 
m m ^ p n -""- 
m r=l 



■■■ ■mm I mini ill llllllliiini mil ■■ n ^Binmin 
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The function of the design algorithm is to solve the last equation by- 
successive iterations- The redesign formula used in each iterative step 

is 

D Ract 

D' = a D - (1-a) — I A u , (B.3.13) 

m m p n r r,m 

'm r=l ' 

where D is the current design variable, and D' represents its improved 
value. The constant a can now be recognized as the relaxation factor , 
which is utilized to control the convergence characteristics of the 
iterative procedure. 

Each iterative design, i.e., each application of (B.3.13) is followed 
by an analysis of resulting structure, so that X and u can be updated. 
The computation of these parameters will be discussed next. 

The displacement gradients u are calculated by the dummy load 
method , which is presented in detail in Reference [2], pp. 37-39. 
Each displacement constraint requires the application of a unit "dummy" 

load to the structure, which must be applied to the same node as the 

* 
displacement constraint u , and in the same direction. If the constraint 

is a rotation, the corresponding dummy load must be unit moment, of 

course. Denoting the element nodal displacement vector due to a dummy 

f rl 
load by {u: }, the gradients of the constrained displacement u are 

u = - I {u. } T [K. ] {u. Cr) } , (B.3.14) 

r.m . i i,m J i ' *■ 

' lem 

where {u. } is the displacement vector of the element caused by the real 

loading . 

Utilizing (B.1.2) and (B.1.3), the derivatives of the element 
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stiffness matrices appearing in (B.3.14) are 



9[K i ] m V 1 m 

[K. ] = -, 1 = n. ([k> J ] + n.A. 1 [k. L J ]) . (B.3.15) 

1 i.m J 3D i v - L x J 11 *■ i JJ *• J 

' m 



Since the unit stiffness matrices of each element are computed and 
stored at the beginning of the program (they do not change upon redesign) , 
the displacement gradients are readily calculated once the solution for 
the dummy loads is available. 

The number of displacement constraints that may be imposed by 
the user is practically unlimited in DESAP 1. Clearly, it would be 
impractical to use a dummy load for each constraint, as this may result 
in a very large dummy load matrix. Moreover, very few of the constraints 
are going to be active in a typical design problem. In view of these 
observations, the maximum number of active displacement constraints is 
assumed to be two in the present version of the program (provision has 
been made to extend this number to four, should the need arise). 

The two candidates for the active constraints are determined after 
the analysis of the current design has been completed; they are chosen 
as the displacements with the two highest displacement ratios 

Q r = V U r * (B.3.16) 

However, if Q < w R , where co < 1 is a pre-determined constant and 
^r max r 

R is the stress ratio defined in (B.2.6), the displacement constraint 
max \ j > f 

is ignored in the next redesign cycle. 

It is important to point out that the derivation of (B.3.14) is 
based on the assumption that the applied loading is size -dependent. 
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Consequently", the formula for the displacement gradients is only a 
rough approximation for thermal stresses and gravity loading. It can 
be shown that the accuracy of the approximation is greater for higher 
values of the inertia exponent n . than for lower values. 

An exact expression for the displacement gradients is available, 
but it requires separate analysis for size-dependent and size- independent 
loading. Since a corresponding method for handling stress constraints 
under size-dependent loads does not exist, it was decided to forgo the 
complications of the exact formula in favor of (B.3.14). 

The Lagrangian multipliers are chosen such as to make the improved 
design displacement-critical, i.e., they are calculated from the con- 
dition u' = u*, r = 1, 2, ...,R .. Using the notation 5u = u' - u 
r r act 6 r r r 

and 6D = D' - D , the change in the nodal displacements can be estimated 
m m m & r 



from the linear approximation 



M 

6u = Z u 6D 
r , r,m m 
m=l 



(B.3.17) 



Passive design variables are governed by the minimum size constraint 

after redesign. Hence, 6D = D* - D , making their contribution to 
6 ' m m m' & 



(B.3.17) 



M 

(5u ) = I u (D*-D ) 
r pass r,m m m 

r m pass 



(B.3.18) 



The improved values of active elements are given by (B.3.13), which re- 
sults in 



(<5u ) + 
r act 



M 
(1-a) Z 
m act 



u D 
r,m m 



R «. 
, act j 

1+— Z X u 
p , s s,m 
m s=l ' i 



(B.3.19) 
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The total change in u is, therefore, 

6u = (6u ) + (Su) . . (B.3.20) 

r r pass r act 

Setting 6u = u* - u , and substituting (B.3.18) and (B.3.19) in 
(B.3.20), we obtain the following simultaneous equations for the 
Lagrangian multipliers X , s = 1,2,..., R ct '- 



R act / M u u ) M 

CI"") £ U s S ^^ D m - -Cl-«) ^ u D n 

s=l \ m act m / m act 



(B.3.21) 



M 

+ E u (D*-D ) - u* + u , r = 1, 2, ..., R . 

r,m m m r r' act 

m pass 

The solution of (B.3.21) requires the prior knowledge of the number 

of active displacement constraints R , and the active-passive identities 

of the design variables. As this information is generally not available, 

the redesign is carried out by an iterative procedure diagramed in 

Fig. B. 3.1. As noted previously, a maximum of two displacement constraints 

are assumed to be active at any stage of the design procedure, i.e., we 

take R + < 2 . 
act - 

The flow diagram inFig.B.3.1 is somewhat simplified, as it omits the 

paths that are followed if the equations for A., and X are singular or 

have a small determinant. In general, if A ■>■ » is indicated by the 

equations, the corresponding displacement constraint is taken as passive, 

i.e., A is treated in the same manner as if it were negative, 
r 
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c 



Compute u 



Start 




Compute D' 
from (B.3.13) 



Label 



D active 

m 



c 



Return 




Yes 



D' => D* . 
m m 

Label 
D passive 



Yes 



Figure B.3.1 
Simplified Flow Diagram for One Redesign Cycle. 
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A design is said to be displacement- critical if 

1 - 6 <Q = max Q <_ 1 + 6 , (B.3.22) 

r 

where Q is the displacement ratio given. in (EL 3. 16), an d 6 is a small 

parameter, also used in defining a stress-critical design see 

(B.2.7). The design is considered acceptable if one of the following 

conditions are met: 

1) All design variables are passive and j< 1 + 6. 

2) The design is displacement-critical and the optimality 
criterion (B.3.12) is satisfied within a prescribed latitude: 



R * 
, act 

55 < — ZXu < - 1 + 5<5 (B.3.23) 

— p , r r,m — v 

K m r= 1 ' 



for all active elements. Experience with the program has shown that the 
"latitude" in (B.3.23) should be considerably less stringent than in 
(B.3.22); hence the use of 56 in the last equation. 

After each analysis, the highest two displacement ratios are printed 
out, and the corresponding load conditions and displacement numbers are 
identified. The quantity 

R . 
, act 

— Z A u , 
p . r r,m 
m r= 1 

which we call the optimality index of the design variable D , is also 

printed for each design variable together with its active-passive 

classification. 



B.4.1 



B.4 Stress and Displacement Constraints 

If both stress and displacement constraints are present, each 
redesign cycle is divided into two parts. First, the structure is 
redesigned with respect to the stresses only; the displacement con- 
straints are ignored in this phase of design. The next step is the 
displacement -cons trained design, where the design variables just obtained 
from the stress-constrained phase are used as the minimum size constraints. 

The displacement-constrained redesign phase will be skipped if the 
analysis of the current design shows that 

< to R , (B.4.1) 

Tiiax — max v -* 

where R and are the maximum stress and displacement ratios de- 
max Tnax c 

fined in (B.2.6) and (B.3.22), respectively. The constant to < 1 is 

also used in conjunction with (B.3.16) in choosing the potentially active 

displacement constraints. 
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B.5 Uniform Scaling Operation 

We recall that the redesign equations for both, stress and dis- 
placement constraints, were based on certain simplifying approximations: 
the stress ratio method assumed that the element nodal forces do not 
change upon redesign; the displacement-constrained redesign formula 
used the linear approximation (B.3.17), and the applied loads were 
assumed to be size dependent. As a consequence of these approximations, 
the redesign process becomes an iterative procedure consisting of re- 
peated applications of the redesign formulas until convergence is 
achieved. 

It is frequently advantageous to interrupt the design procedure by 
the so-called uniform scaling operation, where all the design variables 
are changed by the same scale factor u: 

D' = u D . (B.5.1) 

mm 

The scale factor is calculated from the condition that the scaled design 
{D 1 } should be critical, i.e. 

1 - 6 < max(R ,0 ) < 1 + 6. (B.5. 2) 

— max Tnax — J 

DESAP 1 gives the user the option of using the uniform scaling 
operation whenever the current design is not critical. Once the design 
has been scaled to the critical state (more than one scaling operation 
may be required to achieve this) , the redesign equations will be applied 
in the usual manner. 

This method of design offers two advantages over the use of the 
redesign equations above. Firstly, the use of the scaling operation 
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results in a sequence of critical designs which are very useful in 
monitoring the design process. In particular, the weight comparison of 
the critical designs can be used to terminate the design whenever the 
weight reduction becomes small, or ceases altogether. This is especially 
important when approximate optimization techniques, such as the fully 
stressed design principle, are used. 

The second advantage of using scaling is that it prevents inter- 
mediate designs from departing excessively from the critical state. 
This in turn has a stabilizing influence on the convergence of the 
design process, and in some problems it even makes the difference be- 
tween a convergent process and no convergence at all. 

The scale factor for the stress constraints is simply 

u^ = R , (B.5.3) 

K max v J 

where R is the maximum stress ratio obtained by the stress ratio 
max 

method -- see (B.2.5). The scale factor is exact that is, the 

resulting design will be precisely stress-critical if the internal 

forces remain unchanged upon scaling. This condition is satisfied, 
apart from statically determinate structures, if all element stiffness 
matrices have the form 

[K.] = [k.]A" , (B.5.4) 

where n is common to all elements, and if the loading is size-independent. 

It can easily be shown that if (B.5.4) is satisfied, the dis- 
placement-constrained uniform scaling is also exact, the scale factor 
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being 

P™ - C^) 1 '" , CB.S.S, 

where is the maximum displacement ratio defined in (B.3.22). 

For the problems where the scaling operation is inexact , the 

scale factor for the displacements must be computed from (B.3.17). To 

obtain the scale factor y which would make the r displacement 

critical, we substitute 6u = u* - u and 6D = D' - D = (y -1)D , 

rrr mmmrm 

obtaining 

CD) M 

u* - u = (y l ■'-1) Z u D . 
r r r , r.m m 

m=l 

Solving for the scale factor, we get 

CD) M 

u l J = Cu*-u )/C Z u D ) + 1 . (B.5.6) 

K r r r , r,m m v J 

m=l 

A scale factor is computed for each potentially active displacement 

constraint from (B.5.6), and the maximum value is used for y , i.e. 

\P> = max i/ D) , r = 1,...,R . (B.5.7) 

r 



The factor used in the scaling operation (B.5.1) is the larger of 

(S) . (D) . 
y and y , i.e. 

y = max (y^.y^) . (B.5.8) 



If the structure meets the conditions for exact scaling, the user 
should specify this in the input to DESAP 1, together with the inertia 
exponent n in (B.5.4). The displacement scale factor will then be 
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computed from (B.5.5) and the analysis of the scaled structure will 
be omitted, resulting in a considerable saving in computer time. 
Otherwise, (B.5.6) and (B.5.7) are used to calculate u , and the 
scaled design will be reanalyzed and if not critical, will be scaled 
again. 

DESAP 1 gives the user the option of dispensing with scaling 
altogether. The no-scaling option could be used if the scaling oper- 
ation is inexact and the intermediate critical designs are of no in- 
terest. This could save a substantial amount of computer time, but at 
the expense of losing some control over the convergence of the design 
procedure . 

Scaling should not be used if the loading is dominated by thermal 
or gravity loads and the element stiffness matrices have the form 
[K.] = [k.]A. (linear size-stiffness relationship). It is readily seen 
that under these circumstances a uniform scaling operation is entirely 
ineffective, since it leaves the stresses and displacements unchanged. 

If the uniform scaling option is chosen, the weight of each critical 

design is calculated and the smallest of these weights, W . , is stored. 
6 & mm 

The design procedure is terminated whenever 

(W-W • )/W . > e , (B.5.9) 

v mm' mm — K J 

where W is the weight of the current critical design and e is a small 
prescribed constant. This cut-off criterion prevents further redesign 
operations if these are going to result in a weight increase. 
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The scaling operation is always by-passed if y > 2, in order to 
avoid the corresponding large weight increase. Experience has shown 
that by preceding uniform scaling with a redesign operation, a lower 
weight is obtained than from the scaling-redesign sequence. It should 
be pointed out that a large scale factor would arise only during the 
first design cycle if the initial design is poorly chosen. 
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C. ORGANIZATION OF THE PROGRAM 

C.l Overlay Tree 

The structure of the program is illustrated with the overlay tree 
in Fig. C.l.l Details of the element subroutines are shown separately 
in Figs. C.l. 2 to C.l. 6. 

Additional elements can be added to the program by replacing the 
call NOELEM in ELTYPE with the new element subroutines. The element code 
numbers (MTYPE) 5 and 8 have been reserved for this purpose. 



MAIN 



nzn 



ERROR 



DEVAR 



INPUTJ 



ELTYPE 



nrzi 



UNITWT 



ELMULT 



INPUTD 



INL 



ELSTIF 



MTYPE»1 



TRUSS 



BEAM 



PLANE 



1 f 



S,8 



SHEAR 



SHELL 



BOUND 



See Figs. C.1.2 to C.1.6 



NOELEM 



CLAMP 



CROSS 



SECTOR 



CALBAN 



Element data processing 



■ Executed once per program run ■ 



ERROR 



ADDSTF 



USOL 



MAXD 



PRINTD 



STRESS 



! REARAN 



ELSTFK 



_r 



DPRINT 



MOVED 



DESIGN 



ELTYPE 



MTYPE=1 



TRUSS 



BEAM 



PLANE 



T 



SHEAR 



SHELL 



DPRINT 



5,8 



BOUND 



See Figs. C.1.2 to C.1.6 



NOELEM 



STRSC 



I Element stress recovery and stress ratio redesign 

! 1 — 7= 



DDEF 



MESG 



ARRAN 



FORCES 



X 



DUMAN U 



MOVED 



-s-K- 



Computation of displacemen 
■Executed once with each redesign or scaling cycle 



Figure C.l.l 
Overlay Tree for DESAP 1 
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Figure C.1.2 
Overlay Tree for Truss Element. 



Element Stress Recovery and Stress 

Figure C.1.3 
Overlay Tree for Beam E] 
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Figure C.1.4 
Overlay Tree for Plane Stress Element 
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Figure C. 1.6 
Overlay Tree for Plate/Shell Element. 
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C.2 Element Subroutines 

The element subroutines of DESAP 1 are constructed in such a 
manner as to facilitate modifications by the user. Each element sub- 
routine package consists of two parts: element data processing, and 
stress recovery combined with stress ratio redesign. The element data, 
including the design information, is read in by the. processor sub- 
routines shown by the dashed lines in Figs. C.1.2 to C.1.6. Provision 
has been made for more than one such processor subroutine with each 
element package, the different processors being distinguished by their 
construction code numbers (KODE) . This scheme allows the user to pro- 
vide his own processor subroutines which will best satisfy his special 
design requirements . 

After the data for an element is read in, it is immediately pro- 
cessed and the results placed on auxiliary storage devices. Standard 
computations, such as the formation of the element stiffness matrix, 
transformation of coordinates, etc., are handled by calling the appro- 
priate element computational subroutines , identified by the solid 
lines. No modification by the user is required here, since these sub- 
routines are independent of the construction details and design criteria. 

The processed element data consists of: 

1) Element unit stiffness matrices and load vectors. 

2) Element unit force recovery matrices and force vectors. 

3) Data associated with failure and local instability criteria. 
The first two sets of data are stored in a standard form (common to all 
elements) , whereas the format of the design information is completely 
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flexible it may, for example, contain parameters of the failure 

equations (as in the subroutines presently provided) , or purely 
numerical data in tabular form. 

Following the analysis of the structure, the element nodal forces 
are recovered by subroutine STRSC. The appropriate design subroutine is 
then used to compute the stresses and carry out stress ratio redesign. 
The design subroutines, also identified by the dashed lines in Figs. C.1.2 
to C.1.6, must be compatible with the corresponding processor sub- 
routines that is, if the user supplies his own processor, he must 

also have a matching design subroutine, both labelled with the same 
construction code. 

Note that the boundary element, Fig. C.l.l, is not subject to 
redesign and, therefore, lacks a construction code and design sub- 
routine. 

The user-supplied processor and design subroutines are to be in- 
serted in the place of NOELEM and RETURN, respectively, shown in 
Figs. C.1.2 to C.1.6. 
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C. 3 Storage Requirements 

The common in-core storage blocks used in DESAP 1 are: 

1) COMMON/JUNK/ 300 locations is used for storage of miscellaneous 

data, as its name implies. 

2) COMMON/ELPAR/ 33 locations serves for storage of parameters that 

control the execution of the program. 

3) COMMON/UNITS/ 11 locations contains the assignments (numbers) of 

the input-output units. 

4) COMMON/EM/ 5066 locations is used primarily for processing of 

element data. 

5) COMMON/CONTR/ 23 locations contains control data for the redesign 

operations . 

6) Unlabelled storage area A(n) . This is the main working area of the 
program and its dimension largely controls the allowable size of the 
structure. The capacity of the program can be adjusted by changing 
"n" in the following statements at the beginning of the main program: 

DIMENSION A(n) 

REAL*8 AD (n/2) 

EQUIVALENCE (A(l) , AD(1)) 

MTOT = n . 
The availability of the required storage in A(n) is checked at various 
stages of the program, and if insufficient, an error message is printed 
(subroutine ERROR) and the execution of the program terminated. The 
minimum value of "n" required in each of the subroutines where A(n) is 
used, is listed in Table C.3.3. It is advisable, however, to use the 
largest possible value for "n", since this would reduce the time spent 
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on the transfer of data between the auxiliary storage devices and the 
core. 

7) COMMON/COMPL 798 locations is used only in the plate-shell element 

subroutines . 

The usage of the common blocks in the various subroutines is shown 
in Table C.3.4. 
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Parameter 


Description 


LD 


Number of constrained displacements 




LL 


Number of load conditions 




M 


Max (LL,2) 




MBAND 


Band width of structural stiffness matrix 




NEQ 


Number of degrees of freedom for the structure, i.e., 
number of equations 




NEQB 


Number of equations in a block 
= (n-4*LL)/[(MBAND+M)*4+l] 




ND 


Number of degrees of freedom for an element (see Table 


C.3.2) 


NI 


Dimension of element design information array 
(see Table C.3.2) 




NS 


Number of stresses calculated for an element 
(see Table C.3.2) 




NU 


Number of unit stiffness matrices for an element 
(see Table C.3.2) 




NUMDV 


Number of design variables 




NUMEL 


Number of elements in structure 




NUMNP 


Number of nodal points in structure 




NUMGEO 


Number of different geometric properties 
used for an element type 




NUMTC 


Maximum number of temperatures for which material 
properties of an element type are specified 




NUMFIX 


Number of fixed end-force sets (beam only) 




NV 


Number of unit load vectors for an element 
(see Table C.3.2) 




NW 


Number of initial stress vectors for an element 
(see Table C.3.2) 





Table C.3.1 
List of Parameters that Determine the Storage Requirements, 
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Element 






Parameters 






NU 


NV 


NW ND 


NS 


NI 


Truss 


1 


1 


1 6 


1 


5 


Beam 


2 


2 


2 24 


12 


11 


Plane 


1 


2 


1 12 


15 


10 


Shear 


1 


1 


1 12 


4 


3 


Plate-Shell 


2 


2 


1 24 


6 


5 


Boundary 


1 


1 


1 6 


2 






Table C.3.2 
Values of Storage Parameters for Various Elements. 
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Subroutine 


Minimum Required Value of "n" 


INPUTD 

INL 

ADDSTF 

USOL 

MAXD 

PRINTD 

STRESS 

DERV 
DDESIN 
ELTYPE 


6*NUMNP+12*LL+2*LD 

6*NUMNP+6*LL+2*NEQB*LL 

4* (MBAND+LL) *NEQB+4*LL 

4* (MBAND+M) *NEQB+NEQB 

3*LD+2*NEQB*LL 

6*NUMNP+6*LL+2*NEQB*LL 

3*NUMDV+4*LL+2*NEQB+NEQ 

NUMDV+NEQ+4*NEQB 

6*NUMDV 

NUMDV+10*NUMNP+m, wliere "m" depends on the element type 
(see below) 


Subroutine 


Minimum Required Value of "m" 


TRUSS 

BEAM 

PLANE 

SHEAR 

SHELL 


(2+5*NUMTC) *NUMMAT+2*NUMGE0 
10*NUMGEO+6*NUMMAT+12*NUMFIX 
(2+8*NUMTC) *NUMMAT+5*NUMGE0 
(2+4*NUMTC)*NUMMAT 
(2+7*NUMTC)*NUMMAT 



Table C.3.3 
Core Storage Requirements for A(n) 
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JUNK 


ELPAR 


EM 


UNITS 


CONTR 


COMPL 


MAIN 
CALBAN 
INPUTD 


X 
X 


X 


X 

X 


X 
X 


X 




INL 
ELSTIF 
ADDSTF 


X 




X 
X 








MAXD 

STRESS 

STRSC 


X 
X 
X 


X 


X 








DESIGN 
DDERV 
ARRAN 


X 

X 
X 


X 




X 


X 

X 




FORCES 
DUMAN 
DERV 


X 
X 
X 




X 








DDESIN 
ERROR 
TRUSS 


X 
X 


X 




X 
X 


X 




RUSS 
FTRUSS 
DTRUSS 


X 

X 
X 




X 
X 


X 






BEAM 
TEAM 
NEWBM 


X 
X 

X 


X 


X 
X 


X 
X 






SLAVE 
DBEAM 
PLANE 


X 
X 


X 


X 

X 


X 






PLNAXl 

PLNAX2 

ELAW 


X 
X 

X 




X 
X 


X 
X 






QUAD 
FORMB 
DPLAN1 


X 
X 

X 


X 


X 








DPLAN2 
SHEAR 
PANEL 


X 
X 
X 


X 


X 


X 
X 






FPANEL 

DPANEL 

SHELL 


X 
X 
X 


X 


X 


X 






PLATE1 

QTSHEL 

SLST 


X 
X 




X 
X 


X 




X 

X 


SLCCT 

DSHELl 

BOUND 


X 
X 


X 




X 




X 


CLAMP 


X 




X 


X 







Table C.3.4 
Usage of Common Blocks in Various Subroutines. 
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C . 4 Requ iremen ts for Au xilia ry Storage Devices 

DESAP 1 requires eight sequentially accessible, auxiliary storage 
units (tapes, discs, or drums). The units, together with the required 
storage capacities, are listed below. The storage parameters can be 
found in Table C.3.1. 
II contains permanent storage of the design variable arrays; it also 

functions as a scratch file in the computation of displacement 

gradients 3u /SD . 
6 r m 

No. of locations = 2*NUMDV + Z (2+3*ND) . 

NUMEL 
12 is used for temporary storage of element stiffness matrices (in 

double precision), as a scratch file during the solution of 
simultaneous equations, and for the temporary storage of the 
structural stiffness matrix and dummy load vectors in the calcu- 
lation of displacement gradients. No. of locations=larger of following: 

a) Z [2+ND*(9+2*ND)] , 
NUMEL 

b) |MBAND/NEQB|*(MBAND*NEQB)*2 , 

c) 2*NEQ*LL , 

d) 2*NEQ*(2+MBAND) . 

T3 is used as a scratch file during the solution of simultaneous equations, 

No. of locations = 2*NEQ*(M+MBAND) . 

18 contains permanent storage of displacement number array, element 

design data (including force recovery matrices, thermal stress vectors 
and design information), and the unit weight of each design variable. 
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No. of locations = 6*NUMNP+NUMDV + 14*NELTYP 

+ E [7+ND+(ND*NU+4*NW)*NS*2+NI]. 
NUMEL 

19 used for scratch file in the assembly of the structural stiffness 

matrix, and in the solution of simultaneous equations under the 
real loads . 

No. of locations = the larger of the following: 

a) |(NUMEL) 1/2 [*[ND*(ND*2+9)+2], 

b) |MBAND/NEQB |* (MBAND*NEQB)*2 . 

I 10- -used for temporary storage of the structural stiffness matrix 
and load vectors; also for scratch file during the solution of 
simultaneous equations under the dummy loads. 

No. of locations = the larger of the following: 

a) |MBAND/NEQB|*(MBAND*NEQB)*2, 

b) 4*NEQ. 

I 11- -permanent storage of minimum allowable values of the design 
variables, structural load multipliers, and the actual load 
vectors . 

No. of locations = NUMDV+4*LL+2* (NEQ*LL) . 

I 12 --permanent storage of element unit stiffness matrices, element unit 

load vectors, and the displacement constraints. 

No. of locations = E [7+ND+2*ND* (ND*NU+4NV) ] +NEQ* (1+2*LL) . 

NUMEL 
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The assignments of the auxiliary storage devices, together with 
those of the input-output units (IR = card reader, IW = printer, and 
IP = card punch) > can be changed by altering the assignment numbers at 
the beginning of the Main Program. 
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C.5 Assembly and Solution of Equations 

The assembly of simultaneous equations (structural stiffness matrix 
and load vectors) and their solution, carried out by ADDSTF and USOL, 
respectively, are essentially the same as in the SOLID SAP program. The 
only significant difference is that DESAP 1 uses double precision arith- 
metic, whereas single precision was employed in SOLID SAP. Reconversion 
into single precision is not recommended, even if your machine carries 
above-average number of digits, because numerical errors in the analysis 
tend to be magnified during the redesign cycle. The cumulative buildup 
of error may even lead to a non-convergence of the design process in 
large structures . 

The equilibrium equations are formed in blocks in the unlabelled 
common area A(n) , and stored on auxiliary storage devices. The banded 
structure of the stiffness matrix is exploited throughout the formation 
of the equations and their solution. The size of a block is determined 
automatically by the computer once the bandwidth and the number of load 
vectors have been established. There must be sufficient space in A(n) 
for at least two equations . 

The equilibrium equations are solved by the Gaussian elimination 
procedure, one block at a time. The structural stiffness matrix is not 
destroyed during the solution procedure. 

As can be seen in Fig. C.l.l, the equation solver USOL is used 
twice in each redesign cycle; once for the real loads and once for 
the dummy loads. It may appear at first that considerable computer 
time may be saved by handling the two sets of loads together in a 
single operation, but this is not necessarily true. The problem is 
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that the use of a dummy load for each specified displacement constraint 
(there is no limit on the number of constraints in DESAP 1) may result 
in a very large load matrix. It was considered preferable to limit the 
dummy loads to the potentially active constraints only, although these 
can be determined only after the analysis with real loads has been com- 
puted. Moreover, if no potentially active displacement constraints are 
found, the analysis with dummy loads is skipped altogether. 
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C6 Cutoff Criteria and Restart Option 

The design process is terminated whenever one of the following cutoff 
criteria are satisfied: 

1) The number of design cycles equals a prescribed number NCYCL. 

A redesign cycle is defined as an application of the redesign equations; 
uniform scaling operation, if used, does not count as a redesign cycle. 
This feature enables the user to evaluate the progress of the design 
procedure after a few design cycles and take corrective action if neces- 
sary (particularly, adjust the relaxation parameter a). 

2) The number of successive uniform scaling operations exceeds a 
prescribed number NSCALE. Program termination here indicates that uni- 
form scaling is an ineffective operation and should not be used. 

3) The structural weight of critical design begins to increase, i.e., 

(W - W . )/W . > £ , (C.6.1) 

v mm" mm — 

where the terms are defined in (B.5.9). 

4) The design is acceptable, i.e., convergence is complete. An 
acceptable design satisfies one of the following optimality conditions: 

(i) The design is fully stressed and the displacement con- 
straints are not violated, i.e., 



1 - 6 < R < 1 + S, 1-6 < R . < 1 + S 
— max — — mm — 

and Q < 1 + 6 , 
max — 



(C.6.2) 



where R , R . and Q are the stress and displacement 
max mm Tnax r 

ratios defined in (B.2.6) and (B.3.22j, respectively, and 

6 is a prescribed (small) constant. 
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(ii) The design is displacement-critical, the optimality 

criterion is satisfied for the displacement constraints, 
and stress constraints are not violated, i.e., 



1 - 6 ±Qmax< 1 + 6 > 

R «. 
1 act 

_ i _ 56 < -J- Z A u < - 1 + 56 for all active D , 

— p , r r,m — m 

■m r=l 



(C.6.3) 



and R < 1 + 6 

max — 



At the user's option, DESAP 1 will produce a restart deck just prior 
to the termination of the program. The restart deck contains the last 
values of the design variables together with their minimum allowable 
values. It can be used to replace the original design variable data 
input deck if the user decides to resubmit the program and use the last 
design as the starting point. 
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D. DESCRIPTION OF INPUT DATA 



D.l General Input Data 
I. Heading Card (20A4) 



80 



HED 



HED = Any alphameric statement that is to appear as the first line 
of output . 

II. Program Control Card (415) 



11 



16 20 



NUMNP 


NELTYP 


LL 


NUMDV 



NUMNP = Number of nodal points in the structure. Include special 

nodes used for beam elements (slave nodes and points used to 
define principal axes) and boundary elements (points used to 
define directions of the elements) . 

NELTYP = Number of different element types used in the structure. Count 
each construction code of an element as a separate element 
type. 

LL = Number of separate loading conditions for which the structure 
is to be designed. 

NUMDV = Number of independent design variables in the structure. 

III. Design Control Card (315, 2F10.0, 315, 2F10.0) 



1 


6 


11 


16 


26 


36 


41 


46 


51 


61 70 


NCYCL 


NSCALE 


KSCALE 


DELTA 


EPSIL 


KPUNCH 


KPRINT 


KDISP 


OMEGA 


ALPA 



NCYCL = Maximum allowable number of redesign cycles (see Sec. C.6). If 
NCYCL = 0, analysis of initial design only will be carried 
out. 

NSCALE = Maximum allowable number of successive scaling operations (see 
Sec. C.6). If blank, computer sets NSCALE = 3. 
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KSCALE = Code for uniform scaling operation. 

KSCALE < 0: omit uniform scaling operation. 

KSCALE = 0: carry out uniform scaling whenever the design is 

not critical, and analyse the scaled structure. Used when 

scaling is not exact (see Sec. B.5). 

KSCALE = n > 0, where n is the exponent of A^ in [Ki] = [ki]Aj|. 

Uniform scaling will be carried out whenever the design is not 

critical, but analysis of the scaled structure will be skipped. 

To be used only when scaling is exact (see Sec. B.5). 

DELTA = The small parameter 6 that specifies the latitude of the cut- 
off criteria listed in Sec. C.6. If blank, computer sets 
DELTA = 0.05. 

EPSIL = The small parameter e that specifies the allowable weight in- 
crease of the structure in the weight -cutoff criterion in Sec. C.6. 
If blank, computer sets EPSIL = 0.1. 

KPUNCH = Code for restart card deck. 

KPUNCH j* 0: punch design variable data cards for the last- 
design before program termination occurs (see Sec. C.6). 
KPUNCH = 0: omit punching. 

KPRINT = Code for printout of nodal displacements. 

KPRINT^ 0: print nodal displacements after each analysis of 

the structure. 

KPRINT = 0: omit printout of nodal displacements. 

KDISP = Code for displacement constraints. 

KDISP ^ 0: displacement constraints are being used in the 

current problem. 

KDISP =0: no displacement constraints exist. 

OMEGA = The comparison parameter to < 1 used in selecting the potentially 
active displacement constraints (see Sec. B.4). If blank, 
computer sets OMEGA = 0.8. 

ALPA = The relaxation factor a < 1 used in the displacement-constrained 
" redesign formula (B.3.13). 



Note: the chances that the iterative design procedure con- 
verges uniformly increase with increasing values of ALPA. 
Consequently, it is better to use an ALPA that is too large 
(resulting in under-relaxation) , than one that is too small 
(over-relaxation) . Large values of ALPA, however, may require 
a larger number of design iterations before the cutoff criteria 
are activated. In order to minimize computer time, it is ad- 
visable to start with a "normal" value, which can be estimated 
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from [3] a .= n/(n+l), where n is the inertia exponent that 
dominates the structural stiffness matrix. After a few 
redesign cycles ALPA may be increased or decreased, de- 
pending on the convergence characteristics of the problem. 

IV. Design Variable Data (15, 2F10.0) 

One card for each design variable; cards must be in ascending order of 
design variable numbers. 



16 25 



N 


AOLD 


AMIN 



N 



Design variable number. 



AOLD = Initial value of design variable. If AOLD < AMIN, computer sets 
AOLD = AMIN. 



AMIN 



Minimum allowable value of design variable. 



Automatic data generation if there exists a sequence of design variables 

N = n, n+1, n+2, ... that have identical values of AOLD and AMIN, it is 
sufficient to include only the last card of sequence in the deck. The 
omitted data will be generated automatically. 

V. Nodal Point Data (715, 3F10.0, 15, F10.0) 

One card for each node; cards do not have to be in node number sequence, 
but the last card in the deck must be for the last node. 



11 



16 



21 



26 31 



36 



46 



56 



66 



71 80 
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N = Node number. 

ID = Motion code for nodal displacements (u , u , u ) and rotations 

xyz 
8,0,0) with respect to the global coordinate axes. 

ID(i) = 0: motion is permitted. 

ID(i) = 1: motion is not permitted (the corresponding degrees 

of freedom are eliminated) . 

ID(i) = m, m > 1: node N is connected to a master node m. 

(used for beam elements only see Ch. E) . 



X,Y,Z = Global coordinates of the node. 
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KN = Node number increment used in automatic generation of nodal 
points (see below) . 

KN = 0: automatic generation is not used. 
KN > 0: use automatic generation. 

T = Nodal point temperature. 

Note: temperature distribution in the structure is prescribed 
by nodal temperatures. The reference temperature (i.e. temper- 
ature of the stress-free state) is specified on the element cards. 

Automatic node generation if a series of nodes exists that satisfies 

the following requirements: 

(a) the nodes are equally spaced along a straight line; 

(b) the temperature distribution along the line of nodes is linear; 

(c) the motion codes for all the nodes are identical; 

(d) the node numbers form the sequence n, n+KN, n+2*KN,..., 
then only cards for the first and last node of the series are required. 
The data for the intermediate nodes will be generated automatically. The 
motion code for the series will be taken from the first card , the node 
increment KN from the last card of the series. 

Automatic motion code generation if a particular degree of freedom is 

to be eliminated from a series of nodes, this can be achieved by using 
-1 as the motion code on the first card, and 1 on the last card of the 
series. The computer will automatically use 1 as the motion code for 
the first and all intermediate node cards of the series. 



VI. 



Element Data 



This data differs for various element types and construction codes; it is 
described from Ch. E onwards. 

VII. Structural Load Multipliers (4F10.0) 

One card for each separate load condition; the cards must be in ascending 
order of load condition numbers . 



1 


11 21 


31 40 


A 


B C 


D 


__ 


STP (1) 











STR(l) = Fraction of element load A which is to be added to the load 
condition (structural load multiplier for element load A) . 
STR(2) = Structural load multiplier for element load B. 
STR(3) = Structural load multiplier for element load C. 
STR(4) = Structural load multiplier for element load D. 
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Element loads A, B, C and D are defined with description of element 
data (Ch. E onwards). The element loads represent thermal and gravity- 
loading or distributed loads of prescribed magnitude. The structural 
load multipliers enable the user to choose any fraction of the element 
loads for any load condition. 

VIII. Displacement Constraint Data (214, 12F6.5) 

One card per load condition for each node that has displacement con- 
straints imposed. Cards must be arranged in an ascending order of 
node numbers. The data must end with a blank card . If no displacement 
constraints are used, omit all cards (including the blank card) . 



15 21 27 33 39 45 
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KIJ.S.) 










7- 



Node number. 



= Load condition number. 



Allowable upper limits (absolute values) on nodal displacements 
(u ,u ,u ) and rotations (6 ,6 ,6 ) in the positive global co- 

ordinate directions, and in the negative coordinate directions 



-u 



imposed. 



-u 



y.y 



u 



y 



o 



x 



P5 



y 



:>&^ 



r ) . A blank means that no constraint is 

Note: the right-hand screw 
rule is used to determine 
positive directions for the 
rotations, as shown in Fig. D.l.l. 
The rotations must be specified 
in radians. 



Fig. D.l.l 



Positive Directions of Displacements 
and Rotations . 
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IX. Nodal Load Data (215, 6F10.4) 

One card per load condition for each node that carries concentrated 
loads or moments. Cards must be arranged in an ascending order of node 
numbers. The data must end with a blank card . If no nodal loads are 
present, use the blank card only. 
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N 
L 
R 



Node number. 

Load condition number. 

Concentrated forces (R ,R ,R 1 and moments (M ,M ,M ") in the 

v x' y z J *■ x y z 

global coordinate directions. Positive values denote loads 
in the positive coordinate directions, negative values in the 

negative coordinate directions (see Fig. D.l.l). 
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BAR ELEMENT 



E . 1 General Information 




Fig. E.l.l 
Typical Bar Element. 



A typical bar element is 
shown in Fig. E.l.l. The 
element coordinate axes y and 
z coincide with the principal 
axes of the cross-section, 
whereas x is the centroidal 
axis . 

The axial displacement is 
assumed to be linear in x; 
consequently, the axial force 
P will be constant throughout 
the length of the element. 



Since bar elements do not carry bending or torsion, the nodal rota- 
tions must be suppressed on the Nodal Point Data cards, utilizing the 
motion code ID(6) see Sec. D.l. For the same reason, only zero con- 
centrated moments are permitted on the Nodal Load Data cards. 

The basic element loads consist of gravity loading (in the three 
global coordinate directions), and temperature increase. The temperature 
increase of an element is assumed to be uniform, and is calculated from 



AT = 1/2(T.+T.) - T j. , 
v 1 y ref ' 



(E.l.l) 



where T. and T. are the nodal temperatures (see Nodal Point Data), and 
T . is the reference temperature specified on the element cards. The ele- 
ment load vector due to gravity is obtained by assigning half the 
gravitational force acting on the element to each node. 
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The program can handle temperature -dependent material properties, 
in which case the properties should be listed for two or more tempera- 
tures. The lowest and highest of these temperatures must cover the range of 
temperatures listed on Nodal Point Data cards. The material properties 
of each element are then obtained from the listing by linear interpolation. 

The output from the analysis will consist of the axial force P for 
each element and each load condition. 

The element data deck for each construction code used must start 
with: 
IV A. Element Control Card (615) 
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11 


16 


21 


26 30 


MTYPE 


NUME 


NUMMAT 


NUMGEO 


KODE 


NUMTC 
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NPAR(l) = Code for element type (MTYPE). For bar elements use NPAR(l) = 1. 

NPAR(2) = Number of bar elements with the specified construction code (NUME) 

NPAR(3) = Number of different materials used for the specified construction 
code (NUMMAT) . 

NPAR(4) = Number of different geometric (cross-sectional) properties used 
for the specified construction code (NUMGEO) . 

NPAR(5) = Construction code (KODE) . 

NPAR(6) = Maximum number of temperatures for which material properties are 
given (NUMTC). If blank, computer sets NPAR(6) = 1. 

The remaining element data depends on the construction code used, and is 
described separately below. 
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E.2 Construction Code No. 1 

This construction code is designed mainly for thin-walled bars 
where the wall thickness only is to be changed during the optimization 
process. The size of the element, namely the cross-sectional area A, 
is then related to the moments of inertia by 

I- = j-A, I- = j-A , (E.2.1) 

y J y z z *■ J 

where i - and i - are the unit moments of inertia. 
J y J z 

If P is positive (tensile), the element is redesigned with respect 
to the allowable tensile stress only. The stress ratio redesign formula 
(B.2.3) then becomes 

A'/A = P/P* , (E.2. 2) 

where A' is the improved area and P* = a*A, a* being the tensile strength. 

For negative (compressive) values of P, Johnson's parabolic formula, 
[Sec. CI, Ref. 4] shown in Fig. E.2.1, is used as the failure criterion. 
For slender bars- failure is governed by Euler buckling, in which case 
the stress ratio redesign formula is 



A'/A = -P/P cr , (E.2. 3) 

where P = a A is the Euler buckling load of the current design. For 
cr cr & & 

short bars, the allowable stress is given by Johnsons parabola; the 
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Johnson's parabola: 



a* 
O = CT*f1 £_i 




Euler curve: a ,. = o 

f cr 



L 

r 

a* 
c 



cr 



E 
C 



failure stress 

length of element 

radius of gyration of the cross section 

compressive strength (yield stress or crippling stress) 

2 2 
Cir E/(L/r) = Euler buckling stress 

Young's modulus 

end-fixity coefficient for buckling 

Figure E.2.1 

Johnson's Parabolic Formula 



resulting stress ratio formula can be shown to be 



PP 



A»/A = " p* 



cr 



P*(P -P*/4) 
c*- cr c J 



(E.2.4) 



where P* = a*A and P = a A. The ratio A' /A is computed about both 
c c cr cr r 

principal axes of bending, and the larger of the values is chosen as the 
stress ratio of the element. 

The elements may be guarded against torsional buckling (applicable 
to open sections), or local buckling of the wall by choosing appropriate 
minimum size constraints. 



VI B. Material Properties Data 

A separate data deck, described below, is required for each 
material . The decks do not have to be in a numerical sequence of 
material numbers. 

Material Control Card (215, F10.0) 
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N 



11 20 



N 


NTC 


WT 



= Material number. 



NTC = Number of temperatures for which the properties of this material 
are given. If blank, computer sets NTC = 1. 

WT = Weight-density of the material. 

Material Properties Cards (5F10.0) 

One card for each temperature for which the properties of this material 
are given. Cards must be in ascending order of temperatures. 



1 


11 21 31 


41 50 


T 


| E a | o* 


a* 
c _ 


m _Z nwiTfrl . 1 




i i'in» ^j 


-n 



PMAT(l) = Temperature for which the properties are given (T) . 

PMAT(2) = Young's modulus of elasticity (E) . 

PMAT(3) = Coefficient of linear thermal expansion (a) . 

PMAT(4) = Tensile strength (a*). 

PMAT(5) = Compressive strength (a*). If blank, computer sets PMAT(5)=PMAT(4) 

VI C. Geometric Properties Data (15, 5X, 3F10.0) 

One card for each geometric property; cards do not have to be in 
numerical order. 
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11 


21 31 40 


N 


blank 


AREA 


I- | I- 

y ' z 
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N = Geometric property number. 

AREA = Cross-sectional area of reference section (see below) 
blank, computer sets AREA = 1.0. 



If 



PGE0(1) = Moment of inertia of the reference section about y-axis (I-) 

If blank, computer sets PGE0(1) = 10.0**6, thus eli- y 
minating Euler buckling as a design consideration. 

PGE0(2) = Moment of inertia of the reference section about z-axis (I-) 



If blank, computer sets PGE0(2) = 10.0 
Euler buckling. 



k 6» eliminating 



The Geometric Properties Data is used solely for the computation of 
the unit moments of inertia j- and j-. Therefore, the reference section 
does not have to coincide with the initial design, but can represent any 
acceptable intermediate design. The unit moments of inertia could, of 
course, be read in directly by setting AREA = 1.0 (or blank), in which 
case PGE0 would be interpreted as the unit inertias. 

VI D. Element Load Multipliers (4F10.0) 

Four cards as shown below: 
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Card 1 : x-gravity 

Card 2: y-gravity 

Card 3: z-gravity 

Card 4 : thermal loading 



EMUL = Fractions of the basic element loads (gravity loading in the 
three global coordinate directions and thermal loading) which 
are to be included in the element loads A, B, C and D. 

Element Load Multipliers simply define the element loads A, B, C ~and D. 
Various multiples of these element loads can be added to each structural 
load condition by the use of the Structural Load Multipliers (see Sec. D.l) 
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VI E. Element Data (615, 4F10.0, 15) 

One card for each element; cards must be arranged in an ascending 
order of element numbers . 



1 


6 


11 


16 


21 


26 


31 


41 


51 


61 


71 75 


IEL 


II 


JJ 


IMAT 


I GEO 


IDV 


FRC 


REFT 


ELPYY 


ELPZZ 


INC 



IEL 

II 

JJ 

IMAT 

I GEO 

IDV 

FRC 

REFT 
ELPYY 

ELPZZ 
INC 



= Element number. 

= Number of node i (see Fig. E.l.l). 

= Number of node j . 

= Material number of element. 

= Geometric property number of element. 

=> Design variable number of element. 

= The design variable fraction n. in eqn. (B.1.2) 
computer sets FRC = 1.0. 



If blank, 



= Reference temperature — see Nodal Point Data in Sec. D.l. 

= The end-fixity coefficient C for buckling about y-axis 

see Fig. E.1.2. 

= The end-fixity coefficient C for buckling about z-axis. 

= Node number increment in automatic generation of element 
data (see below). If blank, computer sets INC = 1. 



Automatic element data generation if there exists a series of elements 

IEL = m, m+1, m+2,..., which satisfies the following requirements: 

(a) The node numbers form the sequences 
II = i, i+INC, i+2*INC,. . . 

JJ = j, j+INC, i+2*INC,.. .; 

(b) The rest of the data is identical for all elements of the series, 
then only the last card of the series will be required. The element data 
for the other elements in the series will be generated automatically. 
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E.3 Construction Code No. 2 

This construction code assumes that all the cross-sectional dimen- 
sions are changed by the same proportion upon redesign. Consequently, 
(E.2.1) is replaced by 

I- = j-A 2 , I- = j-A 2 . (E.3.1) 

y J y ' z z *• J 

This scheme allows the user to choose the optimal proportions of the 
cross section (usually determined by local buckling considerations, 
such as crippling stress or torsional buckling) , and assures him that 
these proportions are not changed upon redesign. 

If P is tensile, (E.2.2) remains valid as the redesign formula. 
For compression, the quadratic size-inertia relationships result in the 
following modifications to (E.2.3) and (E.2.4), respectively: 

A'/A = (-P/P cr ) 1/2 , (E.3.2) 

A'/A.= P*/(4P cr ) - P/P* • (E.3. 3) 

The input data and the remainder of the redesign procedure are 
identical to those used in construction code No. 1. 
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BEAM ELEMENT 



F.l General Information 




M- M 
x 



Figure F.l.l 

Typical Beam Element Showing Directions 
of Positive End-Forces and Moments. 



The element coordinate axes y and z shown in Fig. F.l.l, are the 
principal axes of the cross section. The direction of the y-axis is 
defined by the plane of i-j-k, where k is a "third" nodal point that 
must be specified in addition to the nodes i and j . The node k may be 
a conventional nodal point of the structure, or a point used solely for 
defining the principal axes. In the latter case, the motion code on 
the node data card should be used to eliminate all degrees of freedom 
for node k . 

Cubic polynomials in x are used for the lateral displacements, 
whereas the axial displacement is assumed to be linear as in the bar 
element. The bending moments are thus confined to linear functions 
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of x. It follows that the finite element analysis of beam elements is 
exact (apart from numerical errors) for elements that carry concentrated 
nodal forces or moments only. 

The basic element loads are gravity loads in the three global co- 
ordinate directions, and fixed-end forces referred to the element coordinate 
directions. Neither thermal stresses nor temperature -dependent material 
properties are included in the current version of the program. 

Distributed lateral loads are specified in the form of fixed-end 
forces and moments. A fixed-end force (or moment) is simply the end 
force (or moment) acting in a clamped-clamped beam subjected to the 
specified lateral loading. An example on the use of fixed-end forces 
is shown in Fig. 1.2. 
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Figure F.1.2 
Example on the Use of Fixed-End Forces. 
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Fixed-end forces for node j 
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The use of fixed-end forces results in the correct end- forces 
(and moments) in a beam element, but the distribution of the internal 
forces (and moments) is generally inaccurate between the ends. 

Any degree of freedom of an end-node (such as i or j) can be 
connected to the degrees-of-freedom of a "master" node m. If such a 
connection is present, the number of the master node should be used as 
the motion code for the "slave" node (see Nodal Point Data in Sec. D.l). 
Note that this scheme precludes using node number 1 for the master node. 

If the x- displacement of node i is specified as a slave of node m, 
the displacement is expressed in terms of the degrees-of-freedom of the 
master node as follows: 

u . = u + (z.-z )8 - (y.-y )6 . (F.l.l) 

xi xm 1 m ym 1 m zm 

The corresponding formulas for the y and z-displacements can be obtained 
from (F.l.l) by cyclic permutation of x,y,z. 

If the rotation about x-axis of node i is specified to be the slave 
degree -of -freedom, the computer sets 

9 . = 9 . (F.1.2) 

xi xm 

Analagous equations are used for slave rotations about the y and z-axis. 

Slave-master relationship is used to represent various rigid 
connections between two nodes. The edge-reinforced plate in Fig. F.1.3 
can be offered as an example. By specifying all the degrees-of-freedom 
of node j to be the slaves of node m, a rigid link is created between 
the two nodes . The link enforces displacement continuity between the 



F.1.4 



Plate elements 




I 

3 



Rigid link 
Beam element 



Figure F.1.3 
Example o£ Master and Slave Nodes, 



beam and plate element, but at the same time it accounts for the inter- 
element eccentricity e. 

It is important to note that the eccentricity e is kept constant 
throughout the redesign procedure, i.e., no allowance is made for the 
changes in e caused by redesign. 

The printed output following each analysis consists of the end- 
forces and moments shown in Fig. F.l.l. 

The first card of the element data deck for each construction code 
used must be: 
VI A. Element Control Card (615) 
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NPAR(l) = Code for element type (MTYPE). For beam elements use NPAR(l) = 2. 
NPAR(2) = Number of elements with the specified construction code (NUME) . 
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NPAR(3) = Number of different materials used for the specified con- 
struction code (NUMMAT) . 

NPAR(4) = Number of different geometric (cross-sectional) properties 
used for the specified construction code (NUMGEO) . 

NPAR(5) = Construction code (KODE) . 

NPAR(6) = Number of different fixed-end force sets used for the specified 
construction code (NUMFX) . 

The rest of the element data is dependent on the specified construction 
code, and is listed separately below. 
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F.2 Construction Code No. 1 

The size A of the element is taken as the cross-sectional area, and 
the size-stiffness relationships are assumed to be linear: 

T x = Jx A > T y = *9 A > l l = J z A » CF.2.1) 

where I- is the torsional constant of the cross section, I- and I- 
x y z 

are the principal moments of inertia, and j-, j-, j- are the size- 
independent unit values . 

This construction code is intended primarily for the use of thin- 
walled sections where the wall thickness only is to be varied during 
design. Equations (F.2.1) represent a good approximation for closed 
sections, provided that the wall thickness is not too large. For open 

sections, the expressions for I- and I- are still valid, but the 

3 
torsional constant has the form I- = i-A . It must be noted, however, 

xx ' 

that I- is very small in comparison to I- and I- for thin-walled, 
open sections. Therefore, its contribution to the rigidity of the 
structure can be neglected in most cases, i.e. j- = may be used. 

The construction code is also valid for plane bending of beams 
with rectangular cross sections where the width of the cross section 
is subject to design. 

The stress ratio redesign formula is based on the allowable 
stresses only; instability criteria are not used. Local instabilities 
could be guarded against by choosing adequate minimum size constraints. 
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The stresses are calculated at four points at each end section 
of the element, as indicated in Fig. F. 2.1 by points A, B, A 1 and B'. 
Locations of A and B on the cross section are chosen by the user; 
points A 1 and B' are located by the computer to be in "mirror image" 
position of A and B. 
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Circular 
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Figure F.2.1 
Allowable Families of Cross Sections. 
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The input is limited to three families of cross sections: 

(a) Cross sections with at least one axis of symmetry. It is also 
advisable to have the centroid and the shear center coincide, 
since the program takes the center of twist to be the centroid 
of the cross section. The axis of symmetry must be specified 
as the y-axis. 

(b) Zee section with identical flanges. The y-axis must be speci- 
fied as the principal axis that passes through the centroids of 
the flanges as shown in Fig. F.2.1. 

(c) Circular cross section. The orientation of y and z axes is 
completely arbitrary; the computer will automatically rotate the 
axes such that the y-axis will coincide with directic of the 
resultant moment M (see Fig. F.2.1). 

The user must specify the section moduli at points A and B only; the 

computer will automatically generate the moduli for A 1 and B' as shown 

below. 

(a) Symmetric sections: 



(z-) A , - -CZ-) A , (z-) A , - czj) A . (z-) A , = ( z-) A , 



(z-) B , = - C z-) B , cz 5 ) B , - (z z -) B , (z-) B , = (z-) B 



(b) Zee sections: 



(z-) A , = -cz-) A . CZ-) A , = -cz £ ) A . (z-) A , = C z-) A , 



CZy) B , = -CZy) B . CZj) B , = -CZ-) B , (z-) B , = (z-) B . 



(F.2.1) 



(F.2.2) 
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The normal stresses due to bending and axial load are computed 
from the formulas 

a B = + (M-/Z- - M-/Z £ ) , (F.2.3) 

a A = + R-/A , (F.2.4) 

respectively, where the minus sign applies to node i, and the plus sign 
to node j . The shear stress due to torsion is obtained from 

t = M-/Z- . (F.2.5) 

xx *■ J 

The stress ratio redesign formula is based on the modified Von Mises 
yield criterion 

{a/a*) 2 + (t/t*) 2 = 1 , (F.2.6) 

where a = a + a , and a*, x* are the allowable stresses. If a > 0, 

a* is taken as the allowable tensile stress a*; if a < 0, a* represents 

the allowable compressive stress a*. It should be noted that (F.2.6) 

reduces to the original Von Mises criterion if a* = a* and t* = a*/i/5. 
& t c t 

Because only the wall thickness changes during design, the 
section moduli have the form 

Z- = z-A , Z- = z-A , Z- = z-A , (F.2.7) 

x x ' y y ' z z v J 

where z-, z-, and z- are the size-independent unit section moduli. 
Utilizing (F.2.7) and (F.2.3) to (F.2.6) in the failure criterion (F.2.6), 
we obtain for the stress ratio redesign formula the following equation: 



$ - [cP^ * C^) 2 J 1/2 , (F.2.8, 

where a , a , and T are the stresses of the current design A. 

A D 
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VI B. Material Properties Data (15, 5X, 6F10.0) 

One card for each material used; the cards do not have to be in numerical 

sequence. 



1 


6 


11 


21 


31 41 


51 


61 70 


N 


blank 


WT 


E 


1 v 1 a* 


1 a* 
1 c 


1 T* 












■< 


fMAl ID J 







N = Material number. 

WT = Weight-density of the material. 

PMAT(l) = Young's modulus of elasticity (E) . 

PMAT(2) = Poisson's ratio (v) . 

PMAT(3J = Allowable tensile stress (a*). 

PMAT(4) = Allowable compressive stress (a*). If blank, computer sets 
PMAT(4) = PMAT(3). C 

PMAT(5) = Allowable shear stress (t*) . If blank, computer sets 
PMAT(5j = 0.577*PMAT(3). 

VI C. Geometric Properties Data (215, 4F10.0/6F10. 0) 

Two cards for each geometric property; cards do not have to be in 

numerical sequence. 
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31 



41 50 
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KSEC 


AREA 


I- 
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I- 
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I- 
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Card 1 
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n 


21 


31 


41 


51 60 
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z* 
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z* 
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z- B 
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z- B 

z 






PflPO <■"•> 




■5-1 








J ^J 







Card 2 



N 



Geometric property number. 
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KSEC = Code for cross section (see Fig. F.2.1). 
KSEC = 1 or blank: symmetric section. 
KSEC = 2: Zee section. 
KSEC = 3: circular section. 

AREA = Cross-sectional area of reference section. 
If blank, computer sets AREA = 1.0. 

PGE0(1) = Torsional constant of the reference section (I-) . 

PGE0(2) = Moment of inertia of the reference section about y-axis (I-) 

PGE0(3) = Moment of inertia of the reference section about z-axis (I-) 
If KSEC = 3 (circular section), computer sets 
PGE0(3) = PGEO(2). 

PGE0(4) = Section modulus of point A of the reference section for 

torsion (Z^) . 
x 

PGEO(5) = Section modulus of point A of the reference section for 
bending about y-axis (Z-) . 

PGE0(6) = Section modulus of point A of the reference section for 
bending about z-axis (Z-) . 

PGE0(7) = Section modulus of point B of the reference section for 

torsion (Z-) . 
x 

PGE0(8) = Section modulus of point B of the reference section for 
bending about y-axis (Z-) . 

PGEO(9) = Section modulus of point B of the reference section for 
bending about z-axis (Z-) . 

Note: If PGEO(i), i = 4, 5,..., 9 are left blank, the 
corresponding torsional or bending stress is set to zero 
in the analysis. For KSEC = 3 (circular section), stresses 
are calculated at points A and A' only; consequently 
PGE0(6) to PGE0(9) are set to zero by the computer. 
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The Geometric Properties Data is used only for the computation of unit 
inertias and section moduli. The reference section thus does not have 
to be the initial design. 



VI D. Element Load Multipliers (4F10.0) 
Three cards as shown below. 
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21 



31 40 
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B 
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B 


c 


D 






— EMU Li 


,^>. tj 





Card 1 
Card 2 
Card 3 



x- gravity 
y-gravity 
z- gravity 



EMUL = Fractions of basic element loads (gravity loads in the global 
coordinate directions) which are to be included in the element 
loads A, B, C and D. 

Lumped nodal forces are used to represent gravity loads, as was done for 
Bar Elements (see Sec. E.l). Fixed-end forces are not computed. 

The element load multipliers define the contribution of gravity to element 
loads A, B, C and D. Additional contribution to the element loads is 
made by the fixed-end forces; this information is specified separately by 
Fixed-End Force Data and Element Data. 

Specified multiples of element loads A, B, C and D can be added to each 
structural load condition by the use of Structural Load Multipliers 
(see Sec. D.l) . 

VI E. Fixed-End Force Data (15, 6F10.0/5X, 6F10.0) 

Two cards for each set of fixed-end forces used in the analysis; cards 

do not have to be in numerical order. If no fixed-end forces are present 

(NUMFX = 0), omit all cards. 



1 


6 


16 


26 


36 


46 


56 65 


N 


R- 

X 


R- 

y 


R- 
z 


M- 

X 


M- 

y 


M- 
z 


blank 


R- 

X 


R- 

y 


R- 
z 


M- 

X 


M- 

y 


M- 
z 












Jl-l v.. 


'■'•j 




i 



Card 1 : node i 
Card 2 : node j 
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N 



= Fixed-end force set number. 



SFT = Fixed-end forces and moments. Use the sign convention shown 
in Fig. F.l.l. 

VI F. Element Data (715, F10.0, 415, 216, 13) 

One card for each element; cards must be in an ascending order of 

element numbers. 



1 


6 11 16 


21 


26 


31 


36 


46 51 56 61 


66 72 


78 80 


IEL 


I J | K 


IMAT 


I GEO 


IDV 


FRC 


A| B C| D 


NI NJ 


INC 




h* IE(3)-^H 










-«-LC(4)-*H 


h<-JC(12) — H 





IEL = Element number. 

I = Number of node i (see Fig. F.l.l). 

J = Number of node j . 

K = Number of node k. This node defines the direction of y-axis, 
and should not lie on the line of i-j. 

IMAT = Material number of element. 

IGEO = Geometric property number of element. 

IDV = Design variable number of element. 

FRC = The design variable fraction n^ in eqn. (B.1.2). If blank, 
computer sets FRC = 1.0. 

LC(4) = Numbers of fixed-end force sets that are to be included in 
element loads A, B, C and D, respectively (also see Element 
Load Multipliers) . 

NI = End release code for node i. The end release code consists of 
a six digit number, each digit corresponding to an end-force 
as shown below. 



1 



3 4 5 6 



R- 

X 


R- 

y 


R- 

z 


M- 

X 


M- 

y 


M- 
z 



If any of the end-forces is known to 
be zero, due to the presence of a 
hinge or a roller, the digit one must 
be used; otherwise the digit should be 
left blank. 
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NJ = End release code for node j . 

INC = Node number increment in automatic generation of element data 
(see below). If blank, computer sets INC = 1. 

Automatic element data generation if there exists a series of elements 

IEL = m, m+1, m+2,..., which satisfies the following requirements: 

(a) The numbers of the end-nodes form the sequences 

I=i, i+INC, i+2*INC, .... 
J=j, j+INC, j+2*INC,...; 

(b) The rest of the data is identical to all elements of the series, 
including the node number K, 

then only the last card of the series will be required. The element data 
for the other elements of the series will be generated automatically. 



F.3.1 
F.3 Construction Code No. 2 

In this construction code all the dimensions of the cross section 
are scaled by the same factor upon redesign. This scheme gives the 
user an opportunity to choose optimal proportions for the cross section, 
and maintain these proportions throughout the design process. 

The expressions for the bending and torsional properties of the 
cross section now take form 

I- = j-A 2 , I- = j-A 2 , I- = j-A 2 , (F.3.1) 

x J x y y z J z ' *• J 

Z- = z-A 3/2 , Z- = z-A 3/2 , Z- = z-A 3/2 , (F.3. 2) 

xxyy'zz v J 

and the formula for stress ratio redesign becomes 



This equation is solved iteratively for A' /A in the manner described 
below. 

Let r = (A' /A) be the ratio obtained from the previous iteration 
(iteration v) , and r = (A'/A) the improved value predicted by the 
current iteration. If 



2 l a A a B l it \2 

2 >? > CF.3.4) 



the improved value is obtained from the solution of the cubic 

. f v+1,1/2 
equation in (r J : 



H ,1 .2, v,-3 , v+1.3/2 °A, V+K1/2 a B 
/I - Gr*0 (r ) (r ) — r-(r ) - — - = . 



\ a ' 



(F.3. 5) 
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If (F.3.4) is not satisfied, the iteration formula is taken as cubic 

v+1 
equation m r : 



, v+1. 3 



a A v+1 
o^ r 



a B,2 
5*) + 









= 



(F.3.6) 



If t = 0, (F.3.5) yields an exact solution of (F.3.3). 
Similarly, (F.3.6) is an exact redesign equation if a. = 

Pi. 

or a = 0. In both cases only one iterative cycle is needed. Other- 
is 

wise the number of iterations is limited to six in each redesign cycle, 
i.e., the results of the sixth iteration are used for the new design 
even if the iterative solution failed to reach convergence. 

The input data and the remaining details of analysis and re- 
design procedures are identical to those used for Construction Code 
No. 1. 
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G. PLANE STRESS QUADRILATERAL OR TRIANGULAR ELEMENT 
G.l General Information 
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,4(K=L) 




00 



Figure G.l.l 
Quadrilateral and Triangular Plane Stress Elements. 
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The basic plane stress element is the quadrilateral shown in 
Fig. G.l.l(a). The nodes 1 and 2 determine the local Cartesian coordinate 
axes x and y, and the relative orientation of the natural element 
coordinates s and t in the manner indicated on the drawing. Note that 
the natural element coordinates take the values +1 or -1 at the sides of 
the element. All four nodes of the element should be on the same plane. 
If the nodes are not co-planar errors will occur in the computation of 
element properties. The magnitude of the error becomes more serious 
with increasing warp of the element. 

The displacement field within the element is taken in the form 

4 
u(s,t) = Z h i (s,t)u i + hgCs.tJoij + h 6 (s,t)a 2 , 
i=l 



(G.l.l) 



v(s,t) = E h i (s,t)v i + h 5 (s,t)a 3 + h 6 (s,t)« 4 , 
i = l 



where u and v are the components of the displacement field in the 

direction of the local Cartesian coordinate axes. u. and v. denote 

' l l 

the displacement components of node i, and a. are generalized coordinates 
in addition to the nodal displacements. The shape functions in (G.l.l) 
are 

h 1 = I(i-s)Cl-t) h 2 = I(l + s)(l-t) 

h 3 = i(l + s)(l + t) h 4 = I(l-s) (1+t) (G.1.2) 

h = 4^- S > 2 h 6 = T (1 " t)2 • 



G.1.3 
The functions h.. to h. are Lagrangian interpolation polynomials; they 
result in a linear displacement field along each side of the element, 
and thus maintain displacement continuity between the elements . The 
functions h,. and h , on the other hand, represent incompatible displace- 
ment modes, because they cause the sides of the elements to be deformed 
into parabolas . 

According to the examples in [1] , the introduction of the incompatible 
modes can lead to a significant improvement of the analysis for certain 
problems . 

Note that h,. and h, vanish at the nodes. Consequently, a.'s are 
internal degrees of freedom, which can be eliminated at the element level 
by the static condensation procedure. The user of the program, however, 
has the option of suppressing the incompatible modes altogether, in which 
case the computer would set a. = 0, j = 1,2,3,4. 

All the area integrals in the computation of element properties are 
carried out numerically, using Legendre-Gauss quadrature with four 
integration points. It can be shown that the numerical integral is 
exact only for rectangular elements. Since the integration error in- 
creases with the skewness of the element, highly skewed quadrilaterals 
are to be avoided if possible. 

The triangular element used in the program is simply a limiting 
case of the quadrilateral when two of the nodes are allowed to coincide, 
as seen in Fig. G. 1.1(b). It requires no special provisions within 
the program. 

The basic element loads are gravitational forces in the three 
coordinate directions, thermal loading, and in-plane pressure acting on 
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one side of the element (the pressure loading is not allowed for the 
undirectionally reinforced element described under Construction Code 1). 

Temperature -dependent material properties may be used in the pro- 
gram, in which case the material properties should be specified for several 
temperatures. The material properties of each element are obtained 
from the specified properties by linear interpolation. The highest and 
lowest temperatures at which the material properties are prescribed should 
cover the entire range of element temperatures. 

The output of each analysis cycle consists of the in-plane stress 

resultants at the origin of the natural coordinates point in 

Fig. G. 1.1(a), and at points A, B, C, and D. The stress resultants at 
point are printed out twice: once with respect to the local Cartesian 
coordinates x and y, and once with respect to the coordinates x and y 
defined by the angle $ in Fig. G.l.l. The mid-side forces are expressed 
in terms of the directions x., y , etc. shown in Fig. G.l.l (a). The 
convention for positive stress resultants is given in Fig. G.1.2. 




Figure G.1.2 
Positive Stress Resultants 
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By using an appropriate value of the stress -printout option code NS, 
the user can suppress the computation of some or all of the stresses. 

The element data deck for each construction code must start with: 
VI A. Element Control Card 



1 


6 


11 


16 


21 


26 


31 35 


MTYPE 


NUME 


NUMMAT 


NUMTC 


KODE 


n 


NUMGEO 




P^ 




— i\rAK 


{.'J 









NPAR(l) = Code for element type (MTYPE) . For plane stress elements use 
NPAR(l) = 3. 

NPAR(2) = Number of elements with the specified construction code (NUME) 

NPAR(3) = Number of different materials used for the given construction 
code (NUMMAT) . 

NPAR(4) = Maximum number of temperatures for which material properties 
are given (NUMTC). If blank, computer sets NPAR(4) = 1. 

NPAR(5) = Construction code (KODE) . 

NPAR(6) = Code for use of incompatible displacement modes (n) . 
NPAR(6) = 0: use incompatible modes. 
NPAR(6) > 0: suppress incompatible modes. 

NPAR(7) = Number of geometric properties (not used for Construction 
Code No. 2) . 
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G.2 Construction Code No. 1 

This construction code deals with unidirectionally stiffened panels, 
All the cross-sectional dimensions shown in Fig. G.2.1, including the 
stiffener spacing w, are assumed to change by the same proportion during 
redesign. 
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Note: A and I are the cross-sectional area and moment of inertia of 
s s 

the stiffener, respectively. 

Figure G.2.1 
Cross Section of Unidirectionally Stiffened Panel. 



The sheet thickness t is chosen as the size of the element, i.e., it 
is taken as the independent dimension of the cross section. 

For the purposes of analysis, the panel is treated as a homogeneous, 
orthotropic slab with equivalent extensional and shear flexibilities: 
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In (G.2.1), £ and f are the principal material coordinates shown in 
Fig. G.2.2 J N/s, N^ and N^^ represent the average membrane stress 
resultants; and 



t = t + A /w 
s 



is the "average" thickness of the panel 



—j "Effective" panel 



Actual panel 




Figure G.2.2 
Unidirectional ly Stiffened Panel 



The following possible modes of failure are considered in redesign 
operation: 

(1) stresses exceeding their allowable values in the sheet; 

(2) general buckling of the panel; 

(3) buckling of the sheet between the stiff eners; 

(4) failure of stiffeners. 

The stress resultants acting at the middle of the panel (s = 0, t = 0) 
only are used in evaluating these failure modes. 



G.2.3 

Other types of local buckling failure can be handled simply by 
choosing suitable cross-sectional proportions of the panel, i.e., pro- 
portioning the reinforcement such that sheet buckling would always occur 
prior to ( or simultaneously with) the local buckling. The user may even 
optimize the cross-sectional proportions by maximizing the configuration- 
efficiency coefficient of the panel in a manner similar to Ref. [5], p. 107. 

(1) Stress limits 

The modified Von Mises yield criterion 

(oWa*) 2 + (cu/a*) 2 - a.oyca*) 2 + (Wx*) 2 = 1 (G.2.2) 

a y a y Ay 

is used for the stress-constrained design of the sheet, where a* and t* 

are allowable normal and shear stresses, respectively. Different allowable 

stresses may be specified for tension (a*) and compression (a*). Equa-- 

tion (G.2.2) reduces to the original Von Mises yield criterion only if 

a* = a* and x* = a*//T. 
t c 

Using (G.2.2), the stress ratio redesign formula can be shown to be 

t'/t = [(N./N*) 2 + (N^/NS) 2 - N,JW(NXN*) + (N^/N^) 2 ] 1/2 , (G.2.3) 
1 x x' y ' y x y x y J ^ xy xy J J 

where N£ = o"*t, N£ = o*t and N:!u = T*t . 
x y xy 

(2) General buckling of panel 

The panel is treated as a homogeneous, orthotropic plate with simple 
supports. It is assumed that only the compressive stress o^ in the 
direction of the stiff eners influences buckling. 

Since the buckling formula for a general quadrilateral is not 
available, the panel is treated as a rectangle with the "effective" width 
b shown in Fig. G.2.2. The effective width is computed automatically; 
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however, the user has the option of substituting his own value. 

According to Ref. [4], Sec. C.2, p. 26, the compressive buckling 
stress of a stiffened plate is 



x cr 



kir E ,t. 2 



12(l-v 2 ) b 



where the buckling coefficient is given by 



1 + 



„ , EI 

i\'-l s 



N wD 



1 + 



A d 2 I 
s s 



1+0. 88A /(wt) 



]}"., 



In (G.2.5), N is the number of bays between stiffeners , and 



(G.2.4) 



(G.2.5) 



D = 



Et v 



12(1-V ) 
The remaining symbols are defined in Fig. G.2.1. 

For large numbers of stiffeners, we can approximate (N-l)/N 
N = b/w, and obtain for the critical stress resultant 



1 and 



where 



k, = 



(N~) 
x cr 



2(£ 2 ' 



k 1 ir E ,. 
1 t b 



2 4 
12(l-v )b 



(G.2.6) 



1 + 



12(l-v )I 



wt' 



1 + 



A d 2 /I 
. s s 



0.12+0.88t/t/_ 



1/2 



+ 1 



(G.2.7) 



Note that the coefficient k does not change upon redesign since the 
cross-sectional proportions are kept constant. The stress ratio re- 
design formula hence becomes 



t'/t = [-V(N.) cr ] 



1/5 



(G.2.8) 



G.2.5 



(3) Buckling of sheet between stiffeners 
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Figure G.2.3 
Sheet Between Stiffeners. 

The sheet between stiffeners , shown in Fig. G.2.3, is treated as a 
long, simply supported plate of effective width w . The general form 
of the compressive buckling stress for each of the three load conditions 
shown is 



a = k 
cr 



2_ 

TT" E 



12(l-v ) 



/ t_\2 

w 



(G.2.9) 



For simple supports, the buckling coefficients for longitudinal 
compression, transverse compression and shear are, respectively 



k^ = 4.00, k^ = 1.00, k,^ = 5.35. 
x y xy 



(G.2.10) 



If all loads are applied simultaneously, the buckling criterion 
(interaction formula) is taken as 



2 2 
-cW(cu) - cu/(cu) + t^/(t,^) = 1 
x x cr y y cr xy xy cr 



(G.2.11) 
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The interaction formula has an exact theoretical basis only when cu = or 

y 

t^a, = 0; otherwise is an approximation. 

The stress ratio redesign formula obtainable from (G.2.9) to (G.2.11) 
is 



t^_ 

t 



N. - 4(t/t)N. + {[N, + 4(t/t)N.] 2 + [1.495(t/t)N_] 2 } 1/2 

T(WZ) 
x cr 



, (G.2.12) 



where 



(N^) = (cu) t 
v x cr x cr 



Equation (G.2.12) remains valid if N,s and N^ are positive (tensile) 



(4) Failure of stiffeners 
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Figure G.2.4 
Stiffener and Sheet Treated as a Column. 



A stiffener and the attached sheet (see Fig. G.2.4) are treated as a 
column of length "a". The value of "a" is calculated by the computer 
in the manner indicated in Fig. G.2.2, or it may be user-specified. 
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If N/s is positive (tensile), the failure analysis is skipped. For 
compressive loading, the Euler-Johnson failure criterion shown in 
Fig. E.2.1 is used. 

The cross-sectional area of the column in Fig. G.2.4 is 

A = A + wt (G.2.13) 

and its moment of inertia can be shown to be 

I = wt 3 /12 + I + d 2 [wt(A /A) 2 + A (1-A /A) 2 ] . (G.2.14) 

The unit moment of inertia, which is invariant of the value of t, 
can be calculated from 

j = I/t 4 . (G.2.15) 

For long columns that fail by Euler buckling, the critical compressive 
stress resultant is 

N cr = ^=-^L-t 3 , (G.2.16) 

aw a (w/t) 

from which we obtain the redesign formula 

t'/t = C-V N cr )1/3 * (G.2.17) 

Johnson's parabola, which governs the failure of short columns, 
leads to the following quadratic equation in t'/t: 

N*(f/t) 2 + N~(t/f ) -. (N*) 2 /(4N ) = , (G.2.18) 

*~ a. \z c»r 

where N* = a*t, a* being the compressive strength of the stiff ener 
(yield stress or crippling stress) . 
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VI B Material Properties Data 

A separate data deck is required for each material used. The 
decks do not have to be in numerical sequence of material numbers. 
Temperature- Independent Data (215, F10.0) 



1 


6 


11 20 


N 


NTC 


WT 



N 



Material number 



NTC = Number of temperatures for which properties of this material 
are given. If blank, computer sets NTC = 1. 

WT = Weight-density of the material. 

Temperature-Dependent Data (8F10.0) 

One card for each temperature for which the properties of this 

material are given. Cards must be in ascending order of temperatures. 

1 11 21 31 41 51 61 71 80 



T 


1 E 


1 v 


| a | a* 


a* 
c 


T* 


a* 
s 


-<= 






- PMAT(8) 






5»- 



PMAT(l) 
PMAT(2) 
PMAT(3) 
PMAT(4) 
PMAT(5) 
PMAT(6) 
PMAT(7) 
PMAT(8) 



Temperature for which the properties are given (T) 

Young's modulus of elasticity (E) . 

Poisson's ratio (v) . 

Coefficient of linear thermal expansion (a) . 

Tensile strength of sheet (a*) . 



Compressive strength of sheet (a*) . 
Shear strength of sheet (x*) . 



If blank, computer sets 
}► PMAT(6) = PMAT(5), 

PMAT(7) = 0.577 PMAT(5) . 



= Compressive strength (yield or crippling stress) of stiffener 
(a*). If blank, computer sets PMAT(8) = PMAT(6). 
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VI C Geometric Properties Data (15, 6F10.0) 

A separate card is required for each element geometry. The cards 
do not have to be in numerical sequence of geometric property numbers. 



1 


6 




16 


26 


36 


46 


56 65 


N 


TH 


W 


SA 


SI 


D 


WE 



N 

TH 
W 
SA 

SI 



= Geometric property number. 

= Sheet thickness t of the reference cross section (see note below). 

= Distance w between the stiff eners of the reference cross section. 

= Cross-sectional area A of the stiffenerin the reference cross 

s 
section. 

= Moment of inertia I s of the stiffener about its own centroidal 
axis in the reference cross section. 

= Distance d between the centroid of the stiffener and the mid-plane 
of the sheet in the reference cross section. 

= Effective width w e of the sheet between the stiffeners in the 
reference cross section (used in buckling analysis of the sheet). 
If blank, computer sets WE = W. 

The Geometric Properties Data is used only for the computation of 
unit cross -sectional properties, i.e., properties for unit sheet thickness. 
The reference cross section, therefore, can be chosen as any design 
that has the desired cross -sectional proportions; it does not have to 
be the initial design. 
VI D Element Load Multipliers (4F10.0) 

Four cards as indicated below: 



WE 



11 



21 



31 



40 



thermal 


x- gravity 


y-gravity 


z -gravity 


Card 1: 
Card 2: 
Card 3: 
Card 4: 


Element load A 


thermal 


x-gravity 


y- gravity 


z- gravity 


Element load B 


thermal 


x-gravity 


y-gravity 


z- gravity 


Element load C 


thermal 


x-gravity 


y-gravity 


z- gravity 


Element load D 


-*J 


EMUL(4,4) =»- 
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EMUL = Fractions of basic element loads (thermal loading and gravity 
loads in the three global coordinate directions) which are to 
be included in the element loads A,B,C and D. 

Element Load Multipliers simply define the element loads A,B,C and D. 
Various multiples of these element loads can be added to each structural 
load condition by the use of Structural Load Multipliers (see Sec. D.l). 
VI E Element Data (715, 5X, 4F10.0/2F10.0, 215) 

Two cards for each element; cards must be arranged in an ascending 
order of element numbers . 



1 


6 


11 


16 


21 


26 


31 


36 


41 


51 


61 


71 80 


IEL 


I 


J 


K 


L 


IMAT 


IDV 


blank 


FRC 


REFT 


AA 


AB 






IEf* 


n — 


— s»- 

















26 30 



BETA 


EFC 


NS 


INC 



Card 1 



Card 2 



IEL = Element number. 

IE = Node numbers I,J,K and L of the four corner nodes (see Fig. G.l.l) 
Nodes 1 and 2 should define the side that is most parallel 
to the stiffenersas indicated in Fig. G.2.2. For a triangular 
element (their use is not recommended for stiffened panels) , 
the same node number must be used for nodes 3 and 4. 

IMAT = Material number of element. 

IDV = Design variable number of element. 

FRC = The design variable fraction n- in eqn. (B.1.2): A. = n.-D . 
Note that A. , the size of the element, is taken as the 
sheet thickness t. If blank, computer sets FRC = 1.0. 

REFT = Reference temperature (temperature of the stress-free state) 
of the element. 

AA = The effective panel length "a" used in computing the Euler 

buckling load (G.2.16) of the stiffeners. If blank, computer 
calculates "a" in the manner shown in Fig. G.2.2. 

AB = The effective panel width b used in computing the general 

buckling stress (G.2.6) of the panel. If blank, computer calcu- 
lates b in the manner shown in Fig. G.2.2. 
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BETA 



EFC 



NS 



= The angle 3 between the stiffeners and the local x-axis 
(see Fig. G.2.2) . 

= The end-fixity coefficient C used in Euler buckling load 

(G.2.16) of the stiffeners. If blank, computer sets EFC =1.0 
(simple supports) . 

= Stress printout code: 



NS 


Points for which stresses are printed 




(see Fig. G.l.l) 


1 


None 


3 





6 


0,A 


12 


0,A,B,C 


15 


0,A,B,C,D 



If blank, computer sets NS = 3. if NS = 15 and the element is 
triangular, computer sets NS = 12. Note that design of the 
element is always based on stresses at point 0, regardless of 
which printout code is used. 

INC = Node number increment in automatic generation of element data 
(see below). If blank, computer sets INC = 1. 

Automatic element data generation if there exists a series of ele- 
ments IEL = m, m+1, m+2,..., which satisfies the following requirements: 

(a) the node numbers form the sequences 
. I=i, i+INC, i+2*INC,. . ., 

J=j, j+INC, j+2*INC,..., 
K=k, k+INC, k+2*INC,. . ., 
L=l, £+INC, £+2* INC,. . . , 

(b) the remainder of the data is identical for all elements of the 
series , 

then only the last card of the series will be required. The element 

data for the other elements of the series will be generated automatically. 
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G.3 Construction Code No. 2 

Isotropic, homogeneous panels subjected to stress constraints 
only are treated under this construction code. The variable dimension 
of the element (element size) is the panel thickness t. 

Two failure criteria are used in the stress ratio redesign: 

a) Von Mises yield criterion: 

(oya*) 2 + (a 2 /a*) 2 - (a-^a*) (a 2 /a*) = 1 , (G.3.1) 

where a, and a~ are the principal stresses, and a* is the allowable 
normal stress. Different allowable stresses may be specified for tension 
(a*) and compression (a*) . The corresponding stress ratio redesign 
formula is 

t'/t = [(N : /N*) 2 + (N 2 /N*) 2 - (N 1 /N*)(N 2 /N*)] 1/2 , (G.3. 2) 

where N* = a*t, and N , N„ are the principal stress resultants. 

b) Maximum shear stress theory of failure: 

x /t* = 1 , (G.3. 3) 

max 

T being the maximum shear stress and t* its allowable value. This 
max & 

leads to the stress ratio redesign formula 

where (N ) is the maximum shear stress resultant, and N* = x*t. 
s max s 

The value of t'/t is chosen as the larger of (G.3. 2) and (G.3. 4); 
however, the user can suppress the use of (G.3. 4) by setting x* = 
on the material data cards. 



G.3.2 
The redesign formulas are applied to all points of the panel where 

the stress printout is requested (points 0,A,B,C or D shown in Fig. G.l.l) 

This is in contrast to Construction Code No. 1, where point only was 

used for redesign. 

VI B Material Properties Data 

A separate data deck is required for each material used. The decks 

do not have to be in numerical sequence of material numbers. 

Material Control Card (215, F10.0) 



1 



11 21 



N 


NTC 


WT 



N 



= Material number. 



NTC = Number of temperatures for which properties of this material 
are given. If blank, computer sets NTC = 1. 

WT = Weight-density of the material. 

Temperature -Dependent Data (7F10.0) 

One card for each temperature for which the properties of this 

material are given. Cards must be in ascending order of temperatures. 



1 




11 




21 




31 


41 


51 


61 70 


T 


E 




V 


a 


a* 


a* 
c 


T* 














'MAT (7) 


















l 









PMAT(1)= Temperature for which the properties are given (T) 

PMAT(2)= Young's modulus of elasticity (E) . 

PMAT(3)= Poisson's ratio (v) . 

PMAT(4}= Coefficient of linear expansion (a). 

PMAT(5)= Tensile strength (a*). 
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PMAT(6)= Compressive strength (a*). If blank computer sets PMAT(6)=PMAT(5) 

PMAT(7)= Shear strength (t*) . If blank, maximum shear stress theory of 
failure will not be used in redesign. 

VI C Element Load Multipliers (4F10.0) 

Four cards as shown below: 



11 



21 



31 



41 



50 



thermal 


pressure 


x-gravity 


y-gravity 


z-gravity 


Card 1: 
Card 2: 
Card 3: 
Card 4 : 


A 


thermal 


pressure 


x-gravity 


y-gravity 


z-gravity 


B 


thermal 


pressure 


x-gravity 


y-gravity 


z- gravity 


C 


thermal 


pressure 


x-gravity 


y-gravity 


z-gravity 


D 


-e 




EMUL(4,5) 


=9- 





EMUL = Fractions of basic element loads (thermal loading, pressure 
acting on side of element, and gravity loads in the three 
global coordinate directions) which are to be included in the 
element loads A,B,C and D. 

Element Load Multipliers define element loads A,B,C and D. Various 
multiples of these loads can be added to each structural load condition 
by the use of Structural Load Multipliers (see Sec. D.l). 
VI D Element Data (715, F5.0, 3F10.0, 215) 

One car if for each element; cards must be arranged in ascending 
order of element numbers. 



1 


6 


11 16 


21 


26 


31 


36 


41 


51 


61 


71 


76 80 


IEL 


I 


| J K 


L 


I MAT 


IDV 


FRC 


REFT 


PRESS 


BETA 


NS 


INC 






TEC/11 


















Ifc(.4J 







IEL 



Element number. 
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IE = Node numbers I,J,K and L of the four corner nodes (see 

Fig. G.l.l). Nodes 1 and 2 should define the side on which 
pressure load is applied, if any. For a triangular element, 
the same node number must be used for nodes 3 and 4. 

IMA.T = Material number of element. 

IDV = Design variable number of the element. 

FRC = The design variable fraction r\i in eqn. (B.1.2): A^ = HiD m . 
Note that A^, the size of the element, is taken as the panel 
thickness t. If blank, computer sets FRC = 1.0. 

REFT = Reference temperature (temperature of the stress-free state) 
of the element. 



PRESS = Compressive force per unit length (pressure resultant) applied 
to side 1-2 of the element. 

BETA = The angle 3 which defines the x and y axes as shown in Fig'. G.l.l. 
The angle is used only to control the stress printout at point 
(see Sec. G.l) . 

NS = Stress printout code (see Sec. G.2). Note that the design of 

the element is based on the stresses at all the points for which 
stress printout is requested. 

INC = Node increment in automatic generation of element data (see Sec. G.2) 



H. QUADRILATERAL SHEAR PANEL 
H.l General Information 



H.l.l 




Figure H.l.l 
Garvey's Shear Panel Idealization, 



A shear panel, shown in Fig. H.l.l, is assumed to resist only shear 
tractions applied to its edges; it has no rigidity with respect to normal 
edge tractions. The resultant forces of the shear tractions Q. , i = 1,2,3,4, 
must form a self-equilibrating system, so that only one of resultants 
is statically independent. 

All four corner nodes should be on the same plane, otherwise errors 
will arise in the computation of element properties. 

Following NASTRAN [6] , the equivalent nodal forces F and F are 
obtained by "lumping" one half of each adjacent edge force into the corner. 
This method was chosen purely on the basis of convenience, because it 
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results in corner forces that are collinear with the diagonals of the 
element. Again, only one of the corner forces is statically independent, 
the remainder being determined by equilibrium requirements. 

In order to calculate the strain energy (i.e., the stiffness matrix) 
of the element, an assumption must be made regarding the distribution of 
the shear stress t in the panel. We adopted the distribution suggested 
by Garvey [7], which satisfies equilibrium equations, but not compati- 
bility. of the resulting displacement field (except for a parallelogram 
and a rectangle). The assumed shear planes are shown in Fig. H.l.l; the 
magnitude of t is taken to be inversely proportional to the square of the 
distance d from the "baseline" AB. The details of formulating the element 
properties can be found in Ref . [6] . 

Two potential difficulties arise in the use of shear panels: 

a) The state of stress assumed by Garvey can exist only in rectangu- 
lar panel, i.e., only a rectangle can be in a state of pure shear when 
subjected to tangential boundary tractions. It is important, therefore, 
not to use severely skewed shear panels in order to avoid erroneous 
results. 

b) The absence of extensional rigidity will, in general, lead to 
a singular stiffness matrix of the structure, unless each shear panel 
is surrounded by elements that are capable of absorbing the extensional 
corner forces (forces in direction other than F.. and F„) . Bar elements 
are commonly used for this purpose, where the stiffness of each bar 
usually includes the extensional rigidities of the adjacent shear panels. 

Shear panels have traditionally been used in aerospace industry 
to idealize stiffened panels: the shear resistance of the sheet is 



H.1.3 
accounted for by the shear panel, whereas the extentional rigidities 
of the sheet and the stiff eners are lumped into surrounding bar elements . 
This idealization was undoubtedly useful in the days of hand compu- 
tation, but its utility has become dubious with the arrival of the finite 
element method. An orthotropic, plane stress quadrilateral element, such 
as described in Sec. G.2, is not only a more precise representation of 
a stiffened panel, but it also requires less input data and no more computer 
time. 

Shear panel has not been excluded from DESAP 1, partly because it 
may still be a useful element in special applications, and partly in 
recognition of its popularity with aircraft designers. 

Only gravity loads are permitted as the basic element loads. 
Thermal stresses are excluded, of course, since thermal expansion cannot 
be resisted by a state of pure shear. Temperature-dependent material 
properties, however, are allowed. 

The printout that follows each analysis cycle consists of the 
shear flow N. = t . t at each of the four nodes, and the average shear 
flow of the panel , given by 

1 4 

N av = 4 * N i • C"' 1 - 1 ) 

i=l 

The first card of the element data deck for each construction code 

must be: 

VI A Element Control Card (515) 



1 


6 


11 


16 


21 25 


MTYPE 


NUME 


NUMMAT 


NUMTC 


KODE 






NPAR(5)- 




>- 
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NPAR(l) = Code for the element type (MTYPE) . For shear panel use MTYPE = 4. 

NPAR(2) = Number of elements with the specified construction code (NUME) . 

NPAR(3) = Number of different materials used for this construction code 
(NUMMAT) . 

NPAR(4) = Maximum number of temperatures for which the material properties 
are given (NUMTC) . If blank, computer sets NPAR(4) = 1. 

NPAR(5) = Construction code (KODE) . 
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H.2 Construction Code No. 1 

Homogeneous, isotropic shear panels are treated under this 
construction code. The thickness t of the panel is the variable dimension 
(size of the element) . Two design criteria are used in the redesign 
equations : 

(1) The average shear stress T should not exceed its prescribed 
limit T* . The corresponding stress ratio redesign formula is 

t'/t = N /N* , (H.2.1) 

av 

where N is the average shear flow defined in (H.l.l), and N* = x*t. 

(2) The average shear stress should not exceed the buckling stress 
T . For a rectangular panel, the critical shear stress is 

cr 12(l-v z ) I D l 

where b is the unsupported width (shorter dimension) of the panel, and 
k represents the buckling coefficient for shear. The stress ratio re- 
design equation obtainable from (H.2. 2) is 

t'/t = (N /N ) 1/5 , (H.2. 3) 

^ av cr' v J 

where N = x t. 
cr cr 

A good approximation for the buckling coefficient is 

k = C x + C 2 (b/a) 2 , < b/a < 1 , (H.2. 4) 

where a is the longer dimension of the panel, and C. and C are constants 
that depend on the edge support conditions only. The values of these 
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constants for various edge supports are given in Fig. H.2.1. The figure 
also shows the corresponding support codes ISU that are used to specify 
the boundary conditions on element cards. 





a 


C = 5.35 








C 1 = 5.35 
C 2 = 7.25 
ISU = 4 


b 




C 2 = 3.99 
ISU - 1 


1 
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///////;//////// 


fc l =8.99 
£c o =5.72 
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^ISU = 2 








C 1 = 5.35 
C 2 = 5.63 
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C 1 = 8.99 
C 2 = 3.29 
ISU = 3 




////S////////// 


C x = 7.07 
C 2 = 3.91 
ISU = 6 
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Simply supported edge 



//'////'/// Clamped edge 

Figure H.2.1 
Edge Supports for Buckling of Shear Panels. 



Equation (H.2.4) agrees precisely with the theoretical buckling 
coefficients in the limiting cases b/a = and b/a = 1 ; in between these 
values, the approximation is somewhat on the conservative side. 

The unsupported dimensions of the panel, a and b, are user-specified. 
However, if the dimensions are left blank, the computer will use the 
"average" panel dimensions shown in Fig. H.2.2. 
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3(K) 

Actual panel 
Mid-side points 
"Average" panel 



Shear Panel and Its Average Dimensions Used 
in Buckling Analysis. 



VI B Material Properties Data 

A separate data deck is required for each material used. The 

decks do not have to be in numerical sequence of material numbers. 

Material Control Card (215, F10.0) 
1 6 11 20 



N 


NTC 


WT 



N 



= Material number. 



NTC = Number of temperatures for which the properties of this material 
are given. If blank, computer sets NTC = 1. 

WT = Weight-density of the material. 

Temperature-Dependent Data (4F10.0) 

One card for each temperature for which the material properties are 

given. Cards must be in ascending order of temperatures. 
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1 


11 21 


31 40 


T 


E | V 


T* 


H 


- PMAT(4) 


: >- 



PMAT(l) = Temperature for which the material properties are given (T) 
PMAT(2) = Young's modulus of elasticity (E) . 
PMAT(3) = Poisson's ratio (v) . 
PMAT(4) = Shear strength (t*) . 
VI C Element Load Multipliers 

Three cards as shown below. 



1 


11 


21 


31 40 


A 


B 


C 


D 


A 


B 


C 


D 


A 


B 


C 


D 


— ^ 


EMUL 


[4,4)- 





Card 1 
Card 2 
Card 3 



x-gravity 
y-gravity 
z-gravity 



EMUL = Fractions of basic element loads (gravity forces in the three 
global coordinate directions) which are to be included in 
the element load cases A,B,C and D. 

Element Load Multipliers define element loads A,B,C and D. Various 
multiples of these loads can be added to each structural load condition 
by the use of Structural Load Multipliers (see Sec. D.l). 
VI D Element Data ( 15, 3F10.0, 15) 

One card for each element; card must be in ascending order of ele- 
ment numbers . 



1 


6 


11 16 


21 


26 


31 


36 


41 


51 


61 


71 75 


IEL 


I 


J K 


L 


IMAT 


IDV 


ISU 


FRC 


AL 


BL 


INC 






IE (4)— 


>■ 

















IEL 



Element number. 
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I = Node number of node 1 (see Fig. H.2.2). 

J = Node number of node 2 . 

K = Node number of node 3. 

L = Node number of node 4 . 

IMAT = Material number of element. If blank, computer sets IMAT = ] 

IDV = Design variable number of element. 

ISU = Code for edge support conditions used in buckling analysis 
(see Fig. H.2.1). If blank, buckling analysis will be 
suppressed. 

FRC = The design variable fraction H^ in eqn. (B.1.2): A^ = n^Dj,,. 
Note that A^, the size of the element, is taken as the panel 
thickness t. If blank, computer sets FRC = 1.0. 

AL = The longer unsupported dimension "a" of the panel used in 

buckling analysis. If blank, "a" is calculated as shown in 
Fig. H.2.2. 

BL = The shorter unsupported dimension "b" of the panel used in 
buckling analysis. If blank, "b" is calculated as shown in 
Fig. H.2.2. 

INC = Node number increment in automatic generation of element 
data (see below). If blank, computer sets INC = 1. 

Automatic element data generation if there exists a series of 

elements IEL = m, m+1, m+2,... which satisfies the following require- 
ments : 

a) the node numbers form the sequences 

I = i, i+INC, i+2*INC,.. 

J = j, j+INC, j+2*INC,. . 

K = k, k+INC, k+2*INC,. . 

L = I, &+INC, £+2*INC,. . 

b) the remainder of the data is identical for all elements of the 
series, 

then only the last card of the series will be required. The element 

data for the other elements of the series will be generated automatically. 
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I. PLATE QUADRILATERAL OR TRIANGULAR ELEMENT 
I . 1 General Information 



•Mid-side points 



2 (J) 



1(1) 

Figure I. 1.1 

Quadrilateral Plate Element. 



The basic plate-shell element is the plane quadrilateral shown in 
Fig. 1. 1.1. The element is made up of four sub-triangles; the coordinates 
of the internal node (node 5) are obtained by averaging the coordinates 
of the four corner nodes. The local Cartesian coordinates x and y are 
determined by the mid-points of the sides in the manner shown on the 
drawing. If the element is orthotropic, x and y denote the principal 
material coordinates; their direction is specified by the angle (3. 

The membrane behavior of the element is obtained by treating each 
sub-triangle as a constant strain element, and then eliminating the 
degrees -of -freedom of node 5 by matrix condensation. 

The procedure of Clough and Felippa[8] is employed to derive the 
bending properties of the quadrilateral. A so-called LCCT-9 element, 
with nine degrees of freedom, is used for each subtriangle. The 
lateral displacement field of this element is a cubic function, but the 
displacements are partially constrained such that the normal slope along 
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each side is linear. After the four triangles are assembled into the 
quadrilateral element, the degrees of freedom of node 5 are again 
eliminated by static condensation. 

The membrane and bending characteristics of the quadrilateral are 
condensed and stored separately in the form of unit stiffness matrices 
and load vectors as described in Sec. B.l. The procedure of condensing 
the unit stiffness matrices before they are combined into the element 
stiffness matrix is valid only if all the sub-triangles are coplanar. 
Therefore, if the four corner nodes do not lie on the same plane, i.e., 
if the quadrilateral is warped, errors will occur in the computation of 
element properties. 

If a triangle, rather than a quadrilateral, is used in the program, 
a single LCCT-9 element will be used for bending, and a constant strain 
triangle for the membrane action. The triangular element with its local 
coordinate system is shown in Fig. 1.1.2. 



1(1) 




2 (J) 
Mid-side points 
Figure 1.1.2 
Triangular Plate Element. 
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The basic element loads consist of lateral pressure, thermal loading, 
and gravity loads in the three global coordinate directions . The temper- 
ature is assumed to be uniform through the thickness of the plate, i.e., 
thermal stresses caused by temperature gradients in the direction normal 
to the element are neglected. 

Provision has been made for temperature-dependent material properties. 
The properties of each material used should be listed for a range of 
temperatures covering the entire temperature spectrum of the structure . 
Linear interpolation is used by the computer to calculate the properties 
of each element for the specified nodal temperatures. 

The output from analysis contains the stress resultants referred 
to the coordinates x and y. The positive directions of the stress re- 
sultants are shown in Fig. 1.1.3. For a quadrilateral, the stress re- 
sultants are calculated at node 5 (Fig. I. 1.1) by averaging the values 





Figure 1.1.3 
Positive Stress Resultants. 
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for the four sub-triangles. In the case of the triangular element, 
the stress resultants are computed at node 3 in Fig. 1.1.2. 

If the plate elements are used to model a curved surface, the nodal 
rotation about the normal to the surface may be suppressed by the use of 
boundary elements (see Ch. J). 

The element data deck for each construction code used must begin 
with : 
VI A Element Control Card (515) 



1 


6 


11 


16 


21 25 


MTYPE 


NUME 


NUMMAT 


NUMTC 


KODE 


■■»i * 


4PAR(6)- 










—1 



NPAR(l) = Code for element type (MTYPE). For plate elements use 
NPAR(l) = 6. 

NPAR(2) = Number of elements with the specified construction code (NUME) . 

NPAR(3) = Number of different materials used for the specified construction 
code (NUMMAT) . 

NPAR(4) = Maximum number of temperatures for which the properties of any 
one material are given (NUMTC) . 

NPAR(5) = Construction code (KODE) . 

The remainder of element data is listed separately for each construction code. 
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I . 2 Construction Code No. 1 

This construction code treats isotropic, homogeneous plate elements. 
The design variable (size of the element) is the plate thickness t. 

The stress ratio redesign formula is obtained from the modified 
Von Mises yield criterion 

(ojo*) 2 * to Jo*) 2 - o^oJio*) 2 + (Wt*) 2 = 1 , (1.2.1) 

k y Ay Ay 

where a* and t* are the allowable normal and shear stresses referred to 
the x and y axis, respectively. Different values of a* may be used for 
tension (a*) and compression (0*) . 

Note that (1.2.1) represents the original , isotropic Von Mises yield 
criterion if o* = o* and x* = G*//3. In that case, the design would be 
independent of the angle 3 in Figs. 1.1.2 and 1.1.3; consequently, 3 will 
only control the printout of the stress resultants. 

The stress ratio formula obtainable from (1.2.1) is the fourth-order 
equation in t'/t: 

(t'/t) 4 - a(t*/t) 2 + b(t'/t) - c = 0, (1.2.2) 

where the minus sign preceding b is applicable to the upper surface of 
the plate and the plus sign to the lower surface, and 

a = (N./N*) 2 + (N./N*) 2 + (N^/N*) 2 - N^/(N*) 2 

b = 2[(N 5? /N*)(M i? /M*) + (N^/N*) (M^/M*) + (N^/N*) (M^/Mp ] 

(1.2.3) 

- [(VN*)(NU/M*) + (VN*)QWM*)] 
x y y x 

c = (NU/M*) 2 + (NU/M*) 2 ■+ (M^/M*) 2 - (M^/M*) (NU/M*) . 
j\~ y ■^v j 
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In (1.2.3) we used the notation 

N* = to* N* = tx* 

s 

(1.2.4) 
M* = t 2 o*/6 , M* = t 2 T*/6 . 

Equation (1.2.2) is solved by the method of successive iterations, 
the iteration equation being 

, v+1.4 , , , v. , v+1.2 n 
(r ) - (a + b/r ) (r ) - c = , 



i.e., 

v+1 
r 



= [^-(a+b/r V ) + /i(a+b/r V ) 2 + c ] 1/2 , (1.2.5) 



where r is the current value of t'/t, and r represents the improved 
value. The starting value of r is taken as one, and the number of 
iterations is limited to ten. Note that for pure bending action or pure 
membrane state of stress b = 0, in which case, (1.2.5) becomes an exact 
expression for t'/t. 

The redesign equation is applied only to the point where the stress 
resultants are calculated, namely node 5 for a quadrilateral, and node 3 
for a triangle. No local instability criteria are used with this con- 
struction code. 
VI B Material Properties Data 

A separate data deck, described below, is required for each material 
used for this construction code. The decks do not have to be in the 
sequence of material numbers . 
Material Control Card (215, F10.0) 

1 6 11 20 



N 


NTC 


WT 
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N 



= Material number. 



NTC = Number of temperatures for which the properties of this material 
are given. If blank, computer sets NTC = 1. 

WT = Weight -density of the material. 

Material Properties Cards (7F10.0) 

One card for each temperature for which the properties of this 

material are given. Cards must be in ascending order of temperatures. 
1 11 21 31 41 51 61 70 



T 


E 


1 V 


a 


°t 


a* 
c 


T* 


h* 




I 


'MAT (7) 






^- 



PMAT(1)= Temperature for which the material properties are given (T) . 
PMAT(2)= Young's modulus of elasticity (E) . 
PMAT(3)= Poisson's ratio (V) . 

PMAT(4)= Coefficient of linear thermal expansion (a) . 
PMAT(5)= Tensile strength (a*). 

PMAT(6)= Compressive strength (a*). If blank, computer sets PMAT(6)=PMAT(5) 
PMAT(7)= Shear strength (t*) . If blank, computer sets PMAT(7)=0.577*PMAT(5) 
VI C Element Load Multipliers (4F10.0) 
Five cards shown below: 



1 


11 




21 


31 40 


A 


B 


C 


D 


A 


B 


C 


D 


A 


B 


C 


D 


A 


B 


C 


D 


A 


B 


C 


D 






El 




^■! 






'1UL 





Card 1 


lateral pressure 


Card 2 


thermal loading. 


Card 3 


x-gravity. 


Card 4 


y-gravity . 


Card 5 


z -gravity. 
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EMUL = Fractions of basic element loads (lateral pressure, thermal 
loading and gravity loads in the three global coordinate 
directions) which are to be included in element loads A,B,C 
and D. 

Element Load Multipliers simply define the element loads A,B,C and D. 

Various multiples of these element loads can be added to each structural 

load condition by the use of Structural Load Multipliers (see Sec. D.l). 

IV D Element Data (815, 3F10.0) 

One card for each element; cards must be arranged in an ascending 

order of element numbers. 



1 


6 


11 16 


21 


26 


31 


36 


41 


51 


61 


71 80 


IEL 


I 


J K 


L 


IMAT 


INC 


IDV 


PRESS 


REFT 


FRC 


BETA 




- T I"! r A \ w 


















*-^KT t ) 







IEL 
IE 

IMAT 
INC 

IDV 
PRESS 

REFT 

FRC 

BETA 



= Element number. 

= Node numbers I,J,K and L of the four corner nodes (see Fig. 1. 1.1) 
For a triangular element use L = 0. 

= Material number of the element. 

= Node number increment in automatic generation of element data 
(see note below). If blank, computer sets INC = 1. 

= Design variable number of element. 

= Lateral pressure (in the positive z-direction) acting on the 
element. 

= Reference temperature (temperature of the stress-free state) 
of the element. 

= The design variable fraction r).? in (B.1.2): A^ = HiD m - Note 
that the size of the element A^ is the thickness of the plate. 
If blank, computer sets FRC = 1.0. 

= The angle B (in degrees) that defines the x and y axes as 

shown in Figs. I. 1.1 and 1.1.2. Note that the stress resultants 
printed out are referred to the x and y axes. 
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Automatic element data generation if there exists a series of 

elements IEL = m, m+1, m+2,..., which satisfy the following conditions: 

(a) the node numbers form the sequences 

I = i, i+INC, i+2*INC,..., 

J = j, j+INC, j+2*INC,..., 

K = k,-k+INC, k+2*INC,..., 

L = A, A+INC, £+2*INC,..., 

(b) the remainder of the data is identical for all elements of the 
series, 

then only the last card of the series will be required. The element data 

for the other elements of the series will be generated automatically. 



I 
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J. BOUNDARY ELEMENTS 

A boundary element is essentially a line element with specified 
extensional and/or rotational spring constant. The element can be 
used for the following purposes: 

(a) modelling of elastic supports, 

(b) enforcement of specified nodal displacements or rotations in 
a given direction, 

(c) computation of support reactions . 

The direction of the boundary element can be specified in two ways, 
as shown in Fig. J.l. 



Boundary element 





(a) 



GO 



Figure J.l 



Two Options for Specifying Direction n of 
a Boundary Element. 



In method (a) the direction is determined by the structural node N and 
a second node I. The latter may also be a structural node, or a node 
specially created for the purpose. In the latter case, the degrees of 



J. 2 



freedom of I should be suppressed on the node card. In method (b) 
the direction of the element is taken as perpendicular to the lines I-J 
and K-L. Again, the points I,J,K and L may be structural nodes, or 
special points listed on the node cards (with suppressed degrees of 
freedom) . This option is particularly useful when boundary elements 
are used to suppress the normal rotations for thin shells. 

The results of the analysis consist of the axial force and/or 
torque in the boundary element. If the boundary element is an extensional 
spring, the axial force is computed from the formula 



P = kS 



(J.l) 



where k is the spring constant and 6 the computed displacement of node N 

(see Fig. J .2 ) . 

Similarly, the torque in a 
rotational spring obtained 
from 



T = k8 , (J- 2) 

with ' T being the torque and 
6 the computed rotation. 



■J^v. Bound, element 

x - T,e 




Structure 



Figure J. 2 

Positive Displacements and 
Forces Acting on the Structure 



Positive directions of the forces and displacements are defined in Fig. J. 2. 
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If a nonzero displacement 5* is prescribed for node N, the following 
force is applied to the node prior to analysis: 

P* = k6* . (J. 3) 

This applied force will result in the specified displacement only if 
the spring constant k of the boundary element is made much larger than 
the corresponding stiffness of the structure. A specified rotation 0* 
is handled in the analogous fashion. 

A regular elastic support (without a prescribed displacement) is 
obtained by setting 6* = and/or 0* = 0. If, in addition, k is made 
sufficiently large, the boundary element will approximate a rigid 
support. The latter option provides a means of enforcing "skewed" 
boundary conditions due to, for example, inclined roller supports. 

Caution must be exercised in the use of boundary elements. Round- 
off errors in the input data and computation of element properties have 
the effect of introducing springs in directions other than the desired 
orientation. The spring constants of these undesired "elements" are 
proportional to the roundoff errors. Although some of the computational 
errors are minimized by the use of double precision arithmetic, it 
is still important not to make the spring constants of the boundary 
elements too large and to specify the directions of the elements with 
great precision. 
VI A Element Control Card (215) 

1 6 10 



MTYPE 



NUME 



J 



-NPAR(2) H 
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NPAR(l) = Code for element type (MTYPE) . For boundary elements use 
NPAR(l) = 7. 

NPAR(2) = Number of boundary elements (NUME) . 

VI B Element Load Multipliers (4F10.0) 



1 


11 21 


31 40 


A 


B | C 


D 




)-£-■ ■ 


Ul'lULi I'+J 





EMUL = Fractions of specified nodal displacements and rotations which 
are to be included in element load cases A,B,C and D. 

Element Load Multipliers define the element loads A,B,C and D. Various 

multiples of these element loads can be added to each structural load 

condition by the use of Structural Load Multipliers as explained in 

Sec. D.l. For boundary elements, the load multipliers are useful if 

some structural load cases include specified nodal displacements, while 

others do not. 

VI C Element Data (915, 5X, 3F10.0) 

One card for each boundary element; cards must be in ascending order 

of element numbers . 



1 


6 


11 16 21 


26 


31 


36 


41 46 


51 


61 


71 80 


IEL 


N 


I | J K 


L 


KD 


KR 


INC (blank 


SD 


SR 


TRACE 




-^ TV. /T"\ «^_ 


















±L,^J 







IEL 

N 

I 
J 
K 
L 



= Element number. 

= Node number of node N (structural node to which the boundary 
element is attached) . 

= Node numbers I,J,K,L that define the direction of the element. 
If node I only is used as shown in Fig. J. 1(a), set J = K = L = 0. 
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KD = Code for extensional boundary element. 

KD = 0: extensional spring is not used. 
KD = 1: extensional spring is used. 

KR = Code for rotational boundary element. 

KR = 0: rotational spring is not used. 
KR = 1: rotational spring is used. 

INC = Node increment used in automatic generation of element data 
(see below) . 

SD = Specified displacement of node N. 

SR = Specified rotation of node N (radians). 

TRACE = Spring stiffness for extensional and/or rotation spring. If 
blank, computer sets TRACE = 10.0**10. 

Automatic element generation — if there exists a series of elements which 

satisfy the following conditions: 

(a) the element numbers form the sequence 
MEMB = m, m+1, m+2, . . . ; 

(b) the node numbers of successive elements are 
N = n, n+INC, n+2*INC,. 
I = i, i+INC, i+2*INC,. 
J = j, j+INC, j+2*INC,. 
K = k, k+INC, k+2*INC,. 
L = £, JL+INC, £+2*INC, . 



or, J = K = L = 

for all elements if node I only 

is used for direction; 



(c) the remainder of the data is identical for all elements of the 
series; 

then only the last card of the series will be required. The element 

data for other elements of the series will be generated automatically. 
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